ID-A079  683  AIR  FORCE  INST  OF  TECH  RRISHT-PATTERSON  AFB  OH  SCHOO— ETC  F/G  22/3 

APPLICATION  AND  COMPARISON  OF  STABLE  PERIODIC  ORBITS  IN  THE  VIC— ETC<U) 
DEC  79  HD  BEEKMAN 

JNCLASSIFIED  AFIT/GA/AA/79D-2  «  NL 


application  and  comparison  op  stable 

PERIODIC  ORBITS  IN  THE  VICINITY  OF 
LAGKANGIAN  POINT'S  L4  AND  L5  TO 
A  FOUR-BODY  TRUTH  MODEL 


/ 

AFIT/GA/AA/79D-2 


THESIS 

WILLI  AN!  D.  BKEKMAN 
CAPTAIN  USAF 


t 


t-  - 

W 


» 


Approved  For  Public  Release;  Distribution  Unlimited 


AF IT/ GA/A A/ 7  9 0- 2 


APPLICATION  AND  COMPARISON  OF  STABLJ- 
PERIODIC  ORBITS  IN  T! IK  ^ICINITY  OF 

;  L/vGRANGIAN  POINTS  L4  AND  L5  TO 

\  *  * 

\  A  FOUR-BODY  TRUTH  MODEL. 


< /  '■  <  THESIS 

PRESENTED  TO  THE  FACULTY  OF  THE  SCHOOL  OF  ENGINEERING 
OF  THE  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 
AIR  UNIVERSITY  (ATC) 

IN  PARTIAL  FULFILLMENT  OF  THE 
REQUIREMENTS  FOR  THE  DEGREE  OF 
MASTER  OF  SCIENCE 


Pi  William  D.y  Bcckmany7  B.S.,  M.A, 
Captain  USAF 

Graduate  Astronautics 

De  cCBHMK  7  9 

V.  ------ 

(P)  A  7 

Approved  for  public  release;  distribution  unlimited 


Acces 

M.A, 

ntis 

dec:  i 

USAF 

bn* 

Ju.  .1 

| 

E’_- 

0  /  '<=)  c/  ^ 


Contents 


w 


List  of  Figures . . . . 


List  of  Tables . . . . . . . . . ,vii 

List  of  Symbols . . . viii 

Abstract . . . . . . . xi 

I.  Introduction . . . . . . . . . .  1 

Background . . . . .  1 

Problem  and  Scope . . . 11 

II.  Problem  Analysis . 15 

Assumptions, . . . . . . . .  15 

Restricted  Four-Body  Equations  of  Motion..., . .  18 

III.  Wheeler's  Orbit . . . . .  23 

Overview  of  Wheeler '  s  Work . . . . .  23 

Assumptions  and  Coordinate  System .  23 

Constants . ! .  24 

Conversion  to  Truth  Model . . .  26 

Verification  of  the  Truth  Model .  31 

Truth  Model  Application . 36 

IV.  Kolenkiewicz  and  Carpenter's  Orbits . . .  54 

Overview  of  K  &  C 's  Work . . . .  54 

Assumptions  and  Coordinate  System . 54 

Constants . . . . .  56 

Conversion  to  Truth  Model . 58 

Truth  Model  Application.. . 62 

V.  Barkham,  Modi  and  Soudack's  Orbit . 71 

Overview  of  B,  M  &  S's  Work .  71 

Assumptions  and  Coordinate  System........ .  72 

Constants . . .  72 

Conversion  to  Truth  Model.... .  74 

Truth  Model  Application . 75 

VI.  Results . 77 

Discussion . . . . . .  77 

Recommendations . . . 83 

Bibliography . . . . . . . .  89 


ii 


Figure 


Page 


LIST  OF  FIGURES 


1  Initial  Condition  Configuration  for 

Truth  Model . .  ,  ,  . . .  21 

2  Wheeler's  Four- Body  Configuration.., . . . ,.25 

3  Wheeler's  Initial  Condition  Configuration., . . . . . . 28 

4  Transformation  from  Wheeler  Frame  to 

Truth  Frame . 29 

5  Wheeler's  Very  Restricted  Four-Body  1/1 

Periodic  Orbit,  1  Orbit . ,33 

6  Wheeler's  Very  Restricted  Four-Body  1/1 

Periodic  Orbit,  3  Orbits . . . . . .  ,  .  34 

7  Wheeler's  Very  Restricted  Four-Body  1/1 

Periodic  Orbit,  20  Orbits . 35 

8  Wheeler's  Orbit  with  Solar  Eccentricity, 

1  Orbit . 37 


9  Wheeler's  Orbit  with  Solar  Eccentricity, 

3  Orbits . 30 

10  Wheeler's  Orbit  with  Solar  Eccentricity, 

20  Orbits . 39 

11  Wheeler ' s  Orbit  with  Sun-Perturbed  Lunar 

Motion,  1  Orbit . 40 

12  Wheeler's  Orbit  with  Sun-Perturbed  Lunar 

Motion,  3  Orbits . ...41 

13  Wheeler's  Orbit  with  Sun-Perturbed  Lunar 

Motion,  20  Orbits . 42 

14  Truth  Model  of  Wheeler  Orbit,  Calculator- 

Accurate  Lunar  Init.  Cond.  ,  1  Orbit. . 46 

15  Truth  Model  of  Wheeler  Orbit,  Calculator 

Accurate  Lunar  Init.  Cond.,  3  Orbits . ...47 

16  Truth  Model  of  Wheeler  Orbit,  Calculator 

Accurate  Lunar  Init.  Cond.,  20  Orbits., . ,.,,,.48 


EL 


35  No  Drift  Display  of  K  &  C  Orbit  II  Truth 

Model,  3  Orbits,... . . . . . ,..,82 


36  No  Drift  Display  of  B,  MSS  Truth  Model 

Orbit,  1  Orbit . . . . . ,...,,,84 


37  No  Drift  Display  of  M,  M  and  S  Truth  Model 

Orbit,  3  Orbits . . . . . 85 


Figure  Page 

38  No  Drift  Display  of  B,  M  and  S  Truth 

Model  Orbit,  20  Orbits.,...,..,.., . . . .  .,..,..86 

39  Drift  of  Wheeler  Orbit  Over  10  years . . . ,..,.87 

Append ix 

A  Heppenheimer ' s  Four-Body  Configuration., . . . 96 


yi 


LIST  OF  TABLLS 


Coefficients  to  Trigonometric  Serie 
Describing  1<  &  C  Orbits...,.., . 


LIST  or  SYMBOLS 


a^,  a  inertial  coordinate;  frame  axes  of  Earth  -  centered 

truth  model 

a  semimajor  axis  of  the  Sun 

s 

a  angle  of  rotation  of  Earth  -  Moon  barycenter  about  the 

Sun  -  Earth  -  Moon  barycenter 

a  rotation  rate  of  B  about  B 

cm  se:n 


B  Earth  -  Moon  barycenter 

cm 


B  Sun  -  Earth  -  Moon  barycenter 

sere 


e 

s 

f 

s 

G 

L4,  LS 


m 

c 


m 


c 


m 

m 


m 

s 


n 

m 


n 

s 


rotating  coordinate  frame  axes  about  Earth  -  Moon 
barycenter  of  Wheeler  system 

eccentricity  of  the  Sun 
true  anomoly  of  the  Sun 
gravitational  constant 

stable  equilibrium  libration  points  of  the  restricted 
throe-body  problem 

mass  of  the  satellite,  negligible 
mass  of  the  Earth 

mass  of  the  Moon 

mass  of  the  Sun 

m  ,  mass  of  the  Moon 
m 

geocentric  gravitational  constant  for  the  satellite 
geocentric  gravitational  constant  for  the  Moon 
mean  motion  of  the  Moon 
mean  motion  of  the  Sun 


viii 


r 

-  o 


position  vector  in  a  fixed  reference  ellipse 

r  Earth  -  satellite  distance 

c 

r  satellite  --  Moon  distance 

cm 

r  satellite  •-  Sun  distance 

cs 

r  Karth  -  Moon  distance 

m 

r  Moon  -  Sun  distance 

ms 

rp  distance  from  the  Sun  to  the  Earth  -  Moon  bnrvcent 

s 

r  Earth  -  Sun  distance 

s 

t  time 

0  angle  of  rotation  of  Wheeler  axis  about  the  Earth 

Moon  barycenter 

0  rotation  rate  of  the  Moon  about  the  Karth 

v  potential  energy 

to  velocity  vector  in  a  fixed  reference  ellipse 

x^t  y satellite  position  coordinates  of  truth  model 

x  /  y  satellite  velocities  of  truth  model 

c  c 

x  y  satellite  accelerations  of  the  truth  model 

c,  c 

x  ,  y  Moon  position  coordinates  of  truth  model 

m  m 

x  ,  y  Moon  velocities  of  truth  model 

m  J  m 

x  y  Moon  accelerations  of  truth  model 

m,  m 

Xs#  yg  Sun  position  coordinates  of  truth  model 

x  ,  y  Sun  velocities  at  truth  model 

s  s 

..  ~  Sun  accelerations  of  truth  model 

x.  y 

s  s 


XX 


§ ,  n  satellite  posit  ion  coyorcl  i  nates  of  B,  M  &  p  inydel 

§*  satellite  velocities  of  B,  M  &  S  model 


x 


In  this  study,  all  previously  discovered  stable  periodic  orLitf 
about  the  triangular  libration  points  are  tested  on  a  planar 
restricted  four-body  truth  model.  The  truth  model  is  an  algorithm 
developed  from  equations  used  by  T.A .  Hoppenheiim-r  for  colony 
location  used  perturbation  theory  and  the  equation  of  the  center. 
Only  one  of  the  orbits,  developed  by  Wheel  01  using  a  very  re¬ 
stricted  four-body  problem  with  Sun  -  Moon  -  Earth  circular  motion, 
is  found  to  be  relatively  stable  for  at  least  twenty  years.  It  is 
prograde  about  1A  having  a  period  in  resonance  with  the  lunar 
synodic  month.  Two  other  orbits,  one  similar  to  Wheeler’s  and 
one  180°  out  of  phase  found  by  Kol cnhiewicz  and  Carpenter  are 
marginally  stable. 
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APPLICATION  AND  COMPARISON  OF  STABLE 


PERIODIC  ORBITS  IN  THE  VICINITY  OF 
LAGRANGIAN  POINTS  L4  AND  L5  TO 
A  FOUR-BODY  TRUTH  MODEL 
I .  INTRODUCTION 

Background 

Foremost  space  scientists  have  been  looking  in  recent  years 
for  existence  of  stable  periodic  orbits  in  which  to  place  long¬ 
term  satellites  for  weather,  observation,  colonization,  and 
protection.  The  subjects  of  resonance  and  stability  are  closely 
related  to  the  problem  of  evolution  of  the  solar  system.  It  is  a 
physically  involving  problem  and  the  methods  available  to  mathe¬ 
matics  today  seem  unsatisfactory  to  produce  pure  nonlinear  ways 
of  attack.  The  linearization  process  in  both  subjects  is  clearly 
of  doubtful  significance,  so  that,  even  if  very  restrictive, 
numerical  solutions  arc  still  the  best  and  more  valuable  sources 
of  information.  It  is  quite  possible  that  we  know  now  very  little 
more  of  the  entire  problem  that  was  known  to  Poincare',  with  the 
advantage  that  we  can  now  compute  much  faster  and  with  much  more 
precision. 

Hardware  developments  like  the  Space  Shuttle  and  inertial 
Upper  Stage  bring  the  day  closer  when  we  will  be  able  to  use 
these  orbits.  Studies  have  revealed  a  number  of  stable  periodic 
orbits  which  may  exist  to  replace  the  unstable  near-earth  syn¬ 
chronous  orbits  and  give  advantages  such  as  a  continuous 
view  of  the  North  Pole,  which  are  unobtainable  today.  The  area  of 


of  interest  for  these  types  of  orbits  is  cislunar  space,  or  the  area 
of  space  where  satellite  motion  is  affected  by  the  gravitational 
fields  of  the  Sun.  Earth  and  Moon  and  any  other  forces  caused  by 
the  gravitational  attraction  of  Jupiter,  Earth  oblateness,  etc. 
are  negligible.  Many  orbits  have  been  found  in  cislunar  space  but 
their  periodicity  and  stability  in  most  cases  have  been  determined 
using  three-body  or  restricted  four-body  equations  rather  than  a 
truth  model. 

Studies  in  the  area  of  orbital  mechanics  near  L4  and  L5  have 
increased  in  complexity  over  the  years.  Let  and  M^  denote 

three  point  masses  ouch  that  ’Ag  1-1  .  The  masses  move  under  the 
influence  of  their  mutual  gravitational  attractions;  the  force 
between  any  two  masses  is  inversely  proportional  to  the  square 
of  their  distance  and  proportional  to  the  product  of  their  masses. 

It  is  well  known  (Ref  27)  that  there  are  in  this  "three  body 
problem"  five  exact  solutions  in  which  the  three  masses  maintain 
a  constant  configuration  which  revolves  with  constant  angular 
velocity.  An  important  specialization  of  the  three  body  problem 
is  the  restricted  three  body  problem  in  which  is  infinitesimal 
and  M  and  move  in  circular  orbits  around  their  barycenter. 

The  smallness  of  means  that  it  does  not  influence  the  notion 
of  and  For  many  purposes  it  is  convenient  to  describe  the 

motion  of  in  a  coordinate  system  which  is  attached  to  and 
M^.  In  this  rotating  coordinate  system  the  five  Lagrange  solutions 
show  up  as  five  fixed  points  at  which  would  be  stationary  if 
placed  there  with  zero  velocity  (i.e.,  zero  velocity  in  the 
rotating  coordinates) .  If  is  further  known  that,  in  this  rotating 


coordinate  system,  may  describe  small  periodic  orbits  about  the 
Lagrange  solutions.  Glyden  therefore  called  the  points  which 
correspond  to  the  Lagrange  solutions  "centres  of  libration";  they 
are  also  often  referred  to  as  "libration  points"  or  "Lagrange 
points. " 

The  libration  points  are  singular  points  of  the  differential 
equations  of  motion  in  the  restricted  problem  of  three  bodies, 
they  are  also  equilibrium  points  since  the  gravitational  forces 
on  a  mass  placed  in  such  a  point  are  balanced  by  the  centrifugal 
force.  Three  libration  points,  the  col  linear  points ,  arc  found 
on  the  line  connecting  the  two  large  masses;  the  ether  two,  the 
triangular  points,  form  equilateral  triangles  with  the  two  large 
masses.  By  linearizing  the  equations  of  motion  Chari ior  (Ref  9) 
showed  that  there  are  two  classes  of  periodic  infinitesimal  orbits 
around  the  triangular  libration  points,  namely  those  with  short 
period  (period  very  nearly  equal  to  that  of  the  period  of  the  two 
large  masses)  and  those  of  long  period  (the  period  depending  on 
the  mass  ratio  of  the  large  masses) .  Each  of  these  classes  con¬ 
sists  of  concentric,  coaxial  and  similar  ellipses  with  semi -major 
and  minor  axes  in  the  ratio  2:1  for  the  short  period  orbits  and  a 
larger  ratio,  again  depending  on  the  mass  ratio,  for  the  long 
period  orbits.  Plummer  in  1911  (Ref  31)  considered  Charlier's 
problem  in  a  more  general  format  and  from  his  results  some  addi-' 
tional  conclusions  can  be  drawn.  For  a  mass  ratio  of  the  two 
large  masses  smaller  than  1/27,  both  classes  of  orbits  around  the 
triangular  points  can  be  expressed  with  trigonometric  functions; 
these  points  are  therefore  called  stable  libration  points. 


T 


Furthermore,  only  one  of  the  classes  of  orbits  around  the  collincar 
libration  points  can  bo  expressed  in  trigonometric  functions,  the 
other  requi  ring  hyperbolic  functions;  the  col  linear  points  are 
therefore  called  unstable  libr<ation  points . 

The  discovery  in  1906  of  the  first  of  a  group  of  asteroids 
which  appear  to  oscillate  (or,  in  astronomical  terms,  librate) 
around  the  Sun  -  Jupiter  triangular  libration  points,  gave  further 
impetus  to  the  study  of  these  motions.  This  first  discovery  was 
called  Zvchillcs  and  since  subsequent  discoveries  were  also  called 
aftc;r  heroes  from  the  Trojan  group.  Brown  in  1911  (Kef  4) 
considered  the  long  period  orbits  around  the  triangular  libration 
points  by  supposing  finite  amplitudes  of  libration  and  discussed 
in  some  detail  the  dependence  of  period  and  orbit  shape  on 
amplitude.  In  another  paper  (Ref  5)  he  discussed  libration 
orbits  for  mass  ratios  greater  than  1/27.  Willard  in  1913 
(Ref  54)  discussed  the  short  period  orbits,  again  of  finite 
amplitude  and  computed  a  number  of  possible  orbits.  Whereas  all 
this  work  was  based  on  the  restricted  problem  of  three  bodies, 
with  the  discovery  of  more  Trojans  additional  theories  emerged 
which  attempted  to  take  into  consideration  the  actual  physical 
circumstances.  Among  the  first  contributions  were  those  by 
Linders  in  1908  (Ref  28)  and  Smart  in  1918  (Ref  39)  and  in 
1923  (Ref  6)  Brown  published  the  explanation  for  his  theory 
which  was  accurate  enough  to  compute  the  position  of  a  Trojan 
asteroid  within  a  few  seconds  of  arc.  This  theory  was  applied 
numerically  to  Achilles  by  Brouwer  in  1933  (Ref  2)  and  to 
Hector,  which  has  a  particularly  large  libration  amplitude,  by 
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Eckert  in  1933  (Ref  17) .  Since  this  theory  was  numerical  it  had  to  be 
set  up  separately  for  each  asteroid.  A  group  theory  was  outlined  by 
Brown  and  Shook  in  1933  (Ref  8)  in  which  the  interesting  direct  and 
indirect  effects  by  Saturn  were  also  discussed,  Herz  in  1943  (Ref  21) 
carried  out  some  of  the  details  of  Brown  and  Shook's  plan.  Further 
work  concerning  the  motion  of  the  Trojans  was  accomplished  by  VJilkens 
in  1917  (Ref  50) ,  1918  (Ref  51) ,  1926  (Ref  52)  and  1932  (Ref  53) . 

Thuring  in  1930  and  1931  (Ref  46)  considered  again  the 
problem  of  the  long  period  motion,  in  particular  the  dependence 
of  the  period  on  amplitude.  His  subsequent  contributions  in 
1950  (Ref  47)  were  largely  based  on  numerical  work  and  his  1959 
paper  (Ref  48)  was  of  particular  interest  because  of  the  applica¬ 
tion  of  an  electronic  digital  computer.  Thuring 's  claim  of  the 
non-existence  of  long  period  orbits  through  any  orbitrary  point 
was  refuted  by  Rabe  in  1961  (Ref  32)  who  made  a  survey  of  numeri¬ 
cally  computed  long  period  libration  orbits,  expressed  in  Fourier 
Series  expansions.  Rabe  also  discussed  some  aspects  of  the 
stability  of  such  periodic  orbits  and  extended  these  studies  and 
his  survey  in  1962  (Ref  33) ;  similar  work  in  the  same  year  (Ref  34) 
was  done  on  libration  orbits  for  the  Earth  -  Moon  system.  Rabe  has 
developed  the  idea  that  such  periodic  orbits  should  be  used  as 
intermediate  orbits  for  computation  of  real,  nonperiodic  orbits. 
Stumpff  in  1963  (Ref  41)  reconsidered  and  refined  Thuring ' s 
theory,  in  particular  with  respect  to  the  relations  between  long 
period  orbits  with  very  large  amplitudes  around  the  triangular 
libration  points  and  nonperiodic  orbits  in  the  neighborhood  of 
the  collinear  libration  points. 
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The  study  of  libration  points  in  the  Earth  -  Moon  system  was 


initiated  by  Klemperer  and  Benedikt  in  1958  (Ref  25)  .  They  argued 
that  in  analogy  with  the  Trojan  asteroids  there  are  to  be  found  in 
the  combined  gravitational  field  of  the  Earth  and  the  Moon,  two 
areas  in  which  natural  or  artificial  bodies  would  move,  while 
maintaining  a  more  or  le.ss  constant  configuration  with  the  Earth 
and  the  Moon.  Again,  as  was  the  case  with  the  Trojans,  a  natural 
discovery  of  such  a  "cloud"  near  I>4  (the  libration  point  60°  ahead 
of  the  Moon)  was  also  reported.  But  since  then,  the  discovery  has 
been  refuted  although  duct  particles  may  remain  in  the  area  tempo¬ 
rarily  according  to  Rooson  (Ref  35)  and  Wolff  (Ref  56)  in  1907. 

More  recently,  interest,  lias  been  shown  in  the  problem  of  the 
influence  of  the  sun  on  motion  close  to  the  libration  points  of 
the  Earth  -  Moon  system  as  well  as  motion  about  the  Earth  and 
the  Moon  itself.  One  possible  model  for  the  Earth  -  Moon  -  Sun 
system  was  proposed  by  Su-Shu  Huang  in  1960  (Ref  22) ,  who  called 
it  the  "very  restricted  four-body  problem. "  Here  the  Earth  and 
Moon  describe  circular  orbits  relative  to  one  another,  and  their 
center  of  mass  describes  a  circular  orbit  around  the  Sun;  all 
these  orbits  are  Keplerian,  lie  in  a  plane,  and  no  perturbations 
are  considered.  Using  this  model  Huang  studied  the  motion  of  a 
fourth  body  of  an  infinitesimal  mass  in  a  similar  manner  as  in  the 
restricted  four-body  problem.  He  concludes  this  model  gives  a 
general  idea  of  where  the  fourth  body  could  or  could  not  go  under 
given  initial  conditions  when  they  are  no  longer  very  near  the 
Earth.  Columbo  in  1961  (Ref  10)  considered  the  motion  near  L4 
and  L5  under  the  influence  of  the  Sun,  and  the  possibility  of 


stabilizing  it  with  a  solar  sail;  in  another  paper  in  1962  (Ref  11) 
he  gave  a  numerical  analysis  of  the  influence  of  the  Moon  orbit's 
eccentricity. 

Ellis  and  Diana  in  1960  (Ref  10)  on  a  parellel  tack  presented 
some  numerically  computed  libration  orbits  in  the  restricted 
problem.  This  was  extended  by  deVries  and  Pauson  in  1962  (Ref  14) 
by  adding  linearized  equations  of  motion  relative  to  a  stable 
libration  point  in  the  restricted  problem  the  principle  effects  of 
a  fourth  body  representing  the  Sun  as  it  is  related  to  the  Earth 
and  Moon.  Two  linear,  second  order  differential  equations  with 
time  varying  coefficients,  were  obtained  which  could  be  solved  in 
powers  of  the  small  parameter  (mass  of  Sun  divided  by  the  cube  of 
the  Earth  -  Sun  distance) .  The  first  order  solution  and  the  most 
significant  parts  of  the  second  order  solution  wore  obtained  and 
for  a  number  of  different  initial  conditions  this  presented  a 
reasonably  close  agreement  with  numerically  integrated  orbits. 

It  appeared  to  the  authors  that  any  so  called  "stability"  was 
strongly  influenced  by  the  Sun  but  it  also  appeared  possible  to 
choose  the  initial  configuration  of  Earth  -  Moon  -  Sun  and  initial 
conditions  of  the  small  particle  such  that  this  influence  was  small 
enough  for  a  usefully  long  "libration  life"  to  be  possible.  In  a 
subsequent  paper  by  deVries  in  1962  (Ref  15)  the  influence  of  the 
Moon's  eccentricity  was  discussed  and  it  was  found  that,  if  the 
Sun  was  introduced  in  the  consideration  of  motion  near  Earth  - 
Moon  libration  points,  the  Moon's  eccentricity  w’ould  have  to  be 
considered  also.  This  was  to  play  an  important  part  in  the  four- 


body  model  used  in  this  thesis.  Michael  in  1963  (Ref  29)  discussed 


orbit  envelopes  which  depend  on  initial  conditions,  based  on  a 
linearized  analysis  of  the  restricted  problem. 

Using  Huang’s  four-body  model  Cronin  in  1964  (Ref  12)  proved 
that  under  certain  conditions  the  fourth  body  has  a  periodic 
motion,  relative  to  a  rotating  coordinate  frame,  near  each  of  the 
libration  points  of  the  restricted  three-body  problem.  Their 
proof  is  based  upon  assumptions  concerning  the  masses  and  distances 
of  the  bodies  which  are  not  satisfied  by  the  Earth  -  Koon  -  Sun 
system. 

Siferd  in  1965  (Ref  38)  used  Huang's  model  for  the  Earth  - 
Koon  -  Sun  system  to  generate  some  periodic  orbits.  Using  a  nu¬ 
merical  integration  procedure,  the  equations  of  motion  for  the  very 
restricted  four-tody  problem  were  integrated  utilizing  a  digital 
computer  until  some  periodic  orbits  were  obtained.  By  this 
technique  eight  periodic  orbits,  in  the  numerical  sense,  with  a 
respect  to  a  rotating  coordinate  system  were  found.  Three  orbits 
were  around  the  Earth,  three  were  around  the  Moon,  and  two  were 
around  LI.  No  periodic  orbits  near  the  triangular  points  were 
obtained. 

Danby  in  1965  (Ref  13)  investigated  the  influence  of  the  Sun 
on  motion  close  to  the  triangular  points  of  the  Earth  -  Moon  system. 
He  felt  the  very  restricted  four-body  model  inadequate  for  his 
investigation  and  therefore  used  a  model  in  which  the  secular 
perturbations  of  the  moon  due  to  the  Sun  were  retained.  The  results 
may  be  said  to  strengthen  the  hope  that  stable  motion  around  the 
triangular  points  of  the  Earth  -  Moon  system  is  possible.  Other 
investigators  include  Tapley  in  1963  (Ref  43)  and  1965  (Ref  44) 


T" 


who  used  a  mpdel  similar  to  the  very  restricted  four-body  model 

except  the  Moon's  orbit  is  inclined  with  respect  to  the  Earth  -  Sun 

plane.  The  equations  of  motion  for  a  particle  near  the  triangular 

/ 

points  of  the  Vcyf'th  -  Moon  system  are  numercially  integrated  on  a 
digital  computer  for  various  initial  conditions.  One  result 
indicate#’  that  a  particle  placed  initially  at  a  triangular  point 
(L4)#vith  zero  relative  velocity  has  an  envelope  of  motion,  centered 
at  L4,  going  through  a  mode  of  expansion  to  a  value  of  one  Earth  - 
Moon  distance  for  the  major  axis  followed  by  a  mode  of  contraction 
to  a  value  of  1/8  Earth  -  Moon  distance  for  the  major  axis.  The 
envelope  repeats  this  sequence  several  times  during  the  2500-day 
period  investigated.  Feldt  and  Shulrnan  in  1966  (Kef  19)  extended 
the  investigation  of  Tapley  to  5000  days  and  found  that  the 
expansion  -  contraction  of  the  envelope  of  motion  did  not  persist 
due  to  a  lunar  encounter  at  approximately  4000  days.  However, 

Tapley  and  Schutz  in  1968  (Ref  45)  discuss  the  effect  of  the 
constants  used  in  the  model  and  found  that  if  more  accurate  values 
were  used,  an  expansion  -  contraction  motion  persisted  for  over 
8000  days.  Katz  in  1975  (Ref  24)  investigated  numerical  orbits  of  a 
satellite  placed  near  L5  but  all  the  initial  conditions  used  did 
not  have  long-term  stability. 

Wolaver  in  1966  (Ref  55)  used  a  linearized  four-body  approxi¬ 
mation  to  demonstrate  that  proper  use  of  initial  conditions  could 
aid  stability  of  an  orbit  in  the  vicinity  of  L4,  Then  Heppenhcimer 
in  1978  (Ref  20)  developed  realistic  models  incorporating  numerical 
four-body  perturbations  in  a  three-dimensional  analysis.  Although 
he  worked  on  resonant  orbits  about  the  Earth,  his  model  sufficiently 
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mirrored  the  real  world  to  be  used  in  this  thesis  us  the  truth 
model  starting  with  the  proper  initial  conditions  for  stable’ 
periodic  orbits  about  the  triangular  libration  points  discussed  by 
various  authors  using  more  simple  models. 

In  the  search  for  stable  periodic  orbits  about  1.4  and  1,5  for 
space  colony  candidates,  four  were  found.  Schechter  in  1960 
(Ref  36)  concluded  a  stable,  periodic  coplanar  orbit  can  exist 
about  the  Sun  perturbed  Eartli  -  Moon  triangular  point.  The:  model 
used  was  a  three-dimensional  analysis  of  the  long-period  features 
of  four-body  motion  about  L4  where  short-term  period  terras  arc 
removed  from  the  Hamiltonian  via  von  Zicpcl's  method  resulting  in 
a  slowly  varying  Hamiltonian.  He  obtained  an  orbit  with  a  period 
of  28.6  days  with  a  1:1  resonance  with  the  Sun  and  Moon  in  a 
rotating  coordinate  frame.  The  elliptical  motion  is  clockwise 
about  L4  and  has  a  sernimajor  axis  of  approximately  60,000  miles. 
Schechter  demonstrated  that  out-of-plane  motion  is  not  seriously 
excited  by  the  Sun  and  has  a  negligible  effect  on  coplanar  motion. 
It  is  this  coplanar  motion  which  is  the  dominant  factor  as  far  as 
stability  is  concerned. 

But  Kolenkiewicz  and  Carpenter  in  1968  (Ref  26)  confirmed 
Schechter 's  orbit  by  numerically  producing  a  somewhat  larger  orbit 
having  the  same  essential  features  of  the  orbit  by  Schechter.  In 
addition,  a  second  similar  orbit  having  a  phase  difference  of 
180°  was  calculated.  It  is  believed  that  the  discrepancy  in  size 
can  be  accounted  for  by  differences  in  the  models  used  for  lunar 
orbit.  If  Schechter' s  lunar  orbit  was  perturbed  elliptical  instead 


of  circular  an  orbit  on  the  order  of  50%  larger  would  result  from 


his  equations. 


Wheel ei  in  1973  (Rof  49)  found  a  stable  periodic  orbit  about. 

L4  in  the  restricted  problem  of  four  bodies.  It  exhibited  a  1:1 

'f- 

coiamonsurabi  lily  with  tlie  Moon's  synodic  month  and  has  a  per iod 

of  29. 9303332  days.  The  eiyenva lues  were  solved  for  and  the 

* 

Poincare'  exponents  were  then  determi nod  to  be  pure  imaginary. 

This  Anplied  linear  stability  of  the  orbit  in  the-  restricted 
four-body  problem. 

Barkham,  Modi,  and  Soudack  in  1975  (Kef  1)  found  a  2;] 
literal  solution  to  the  restricted  four-body  problem  about  l4  and 
L5  and  numerically  generated  periodic,  four-body  solutions  that 
agreed  to  within  5t  of  the  literal  solution. 

Last  of  all,  Kamel  and  Breakwell  in  1970  (Kef  23)  found 
similar  results  to  Kolcnkiowi.cz  and  Carpenter  using  the  von  Zei.pel 
technique.  More  surveys  of  motion  near  the  triangular  Earth  - 
Moon  libration  points  are  gi.ven  by  deVries  (Ref  16)  ,  Stag  and 
deVries  (Rof  40) ,  Szebehely  (Ref  42) ,  and  Schutz  (Ref  37) . 

Problem  and  Scope 

With  space  colonization  on  the  horizon,  it  appears  that  the 
most  fundamental  questions  about  motion  near  libration  points  are 
those  about  the  existence  of  periodic  orbits  and  the  stability  of 
such  orbits.  If  stable  periodic  solutions  exist,  solutions  of 
differential  equations  at  or  near  conditions  of  common surabil ity 
may  be  used  as  intermediate  orbits  for  the  computation  of  non¬ 
periodic  orbits  by  perturbation  analysis.  In  the  restricted 
problem  of  three  bodies  the  existence  of  periodic  orbits  about,  the 
triangular  libration  points  is  well  established.  This  result 


of  motion 


followed  from  the  analysis  of  the  linearized  equation;.; 
and  served  to  exhibit  the  stability  of  the  triangular  configuration 
as  one  of  Lagrange's  exact  solutions  of  the  restricted  problem,  onl 
in  so  far  as  the  linearization  is  valid,  that  in,  only  for  infini- 
testimal  disturbances.  The  apparent  existence  of  non- infinitesimal 
periodic  orbits  (Brown,  Thuring,  Rabe)  followed  either  from  the 
analysis  of  higher  order  approximations  of  the  differential 
equations  (but  still  not  exact)  or  from  numerical  work.  It  is 
very  difficult  to  derive  meaningful  results  by  qualitative  methods 
and  with  the  problem  of  libration  orbits  we  may  still  be  in  the 
position  of  trying  to  come  to  specific  results  by  the  study  of 
particular  analytical  or  numerical  solutions. 

Considering  the  modern  trends  in  the  study  of  nonlinear 
mechanics  toward  qualitative  methods  one  may  expect  that  any  new 
work  on  triangular  libration  points  should  concentrate  on  the 
establishment  of  a  proof  of  stability  of  libration  orbits.  If 
then  a  solution  in  the  form  of  analytical  expressions  pf  the 
coordinates  as  functions  of  time  with  an  exhibition  of  integration 
constants  would  be  at  all  required,  one  should  use  periodic 
orbits  (whose  existence  would  first  bo  proved)  as  intermediate 
orbits  for  the  perturbation  analysis.  Two  reasons  discourage 
one  from  the  following  approach.  First  of  all,  even  though  the 
past  fev;  decades  have  seen  a  significant  development  of  methods 
and  theorems  in  nonlinear  mechanics  there  is  still  very  little 
known  about  systems  of  higher  than  second  order.  The  methods 
of  the  phase  plane,  so  convenient  and  easily  visualized  for 
second  order  systems,  must  be  transferred  to  multidimensional 
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phase  space  which  introduces  some  formidable  complications. 
Secondly,  the  few  qualitative  results  which  are  known  about  the 
triangular  libration  points  specifically  have  been  derived  only 
for  the  restricted  three-body  problem  which  is  really  very  special 
since  its  Hamiltonian  docs  not.  contain  the  independent  variable 
explicitly.  On  the  other  hand,  preliminary  studies  have  shown 
clearly  that  -in  the  use  of  the  1,'arth  -  Moon  libration  points  the 
influences  of  the  Sun  as  the  fourth  bodv  and  of  the  Moon's  orbital 
eccentricity  are  quite  important.  The  Hamiltonian  of  such  a  problc 
contains  the  independent  variable  in  periodic  terms  of  short  and 
long  periods,  and  especially  with  periods  commensurable,  or  nearly 
so,  with  the  principle,  periods  of  the  problem.  Very  little  is 
known  at  all  about  how  certain  qualitative  results  derived  for 
constant  Hamiltonian  could  be  transferred  to  a  similar  problem 
with  time  varying  Hamiltonian. 

So  finding  a  solution  to  the  full  analytical  equations  of 
motion  describing  real  world  forces  is  not.  possible  to  date.  The 
best  method  to  treat  the  problem  of  colony  location  is  to  first 
develop  reference  colony  orbits  in  the  restricted  three-body 
problem.  Such  orbits  arc  then  studies  in  a  very  restricted  four- 
body  problem,  wherein  the  motions  of  the  Earth,  Moon,  and  colony 
are  determined  by  numerical  integration.  This  author  has  found 
four  orbits  around  L4  and  1,5  claimed  to  be  stable  and  periodic 
using  very  restricted  four-body  equations.  This  thesis  will 
attempt  to  test  those  orbits  using  a  model  more  closely  resembling 
the  real  world  to  see  if  the  orbits  are  stable,  and  similar.  If 


they  are  not  stable,  this  paper  will  attempt  to  explain  why  they 


arc  not. 


The  Layi  an  jiai;  point  s  J,4  and  1,5  will  be  used  as  reference 
points.  They  are  the  stable  equilibrium  points  in  the  restricted 
tlirec—  body  problem  for  mass  ratios  less  than  .0385.  Here  the 
satellite  remains  fixed  relative  to  the  other  two  bodies,  if  given 
the  correct  initial  velocity.  However,  in  cislunar  space  L4  and  L5 
are  no  longer  equilibrium  points  although  we  will  still  refer  to 
colony  motion  about  these  points.  Any  truth  model  used  should 
reduce  to  the  restricted  three-body  equations  of  motion.  Li  and 
L5  passes  a  triangular  symmetry  with  each  other  in  relation  to  the 
Earth  and  Moon  so  that  motion  about  1,4  can  be  considered  identical 
to  motion  about  L5.  Exceptions  to  this  are  when  perturbative 
effects  of  other  planetary  bodies  are  taken  into  account.  Then 
the  motion  about  L4  will  exhibit  slightly  different  solutions 
compared  to  motion  about  L5.  This  author  neglects  tho.no  minute- 
planetary  perturbations  so  when  L4  is  referred  to,  it  can  be 
considered  L5  as  well. 
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1 1 •  PROBLEM  analysis 

Assumptions 

The  stable  periodic  orbits  about  L4  found  to  date  were  computed 
using  very  restricted  four -body  (VRFB)  models.  The  VRFB  models 
neglect  the  important  indirect  effect  of  the  Sun,  i.e.,  the 
gravitational  effect  on  the  motion  of  the  Earth  and  Moon,  Hence, 
the  results  obtained  from  the  simplified  VRFB  model  cannot  be  used 
to  infer  motion  in  the  real  world.  That  is,  based  on  VRFB  results, 
it  is  not  known  whether  the  expansion  -  contraction  of  the  envelope 
of  motion  which  exists  in  the  VRFB  models  also  exists  in  the  real 
world.  Instead  the  nature  of  the  solar  influenced  particle  motion 
near  L4  will  be  studied  by  numercially  integrating  the  equations 
of  motion  over  a  number  of  years  using  a  model  which  closely 
represents  the  real  world. 

The  model  used  to  test  these  orbits  is  one  used  by 
T. A.  Heppenheimer  (Ref  20)  in  a  paper  locating  space  colonies 
in  high  Earth  orbits.  In  using  this  model  the  assumptions  made 
are: 

1)  The  Sun,  Earth  and  Moon  are  considered  to  be  point  masses. 

2)  The  mass  of  the  colony  satellite  is  negligible  compared 
to  the  masses  of  the  other  three  bodies  and  exerts  no  forces  to 
affect  their  motion. 

3)  The  gravitational  forces  of  other  planets  have  a  negligible 
effect  on  motion  and  is  ignored. 

4)  The  motion  of  all  four  bodies  is  limited  to  one  plane. 

5)  The  motion  of  the  Sun  is  taken  to  be  an  unperturbed  ellipse 
with  respect  to  the  Earth  -  Moon  barycenter. 
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The  satellite  motion  about  L4  is  far  enough  away  from  the 
three  other  attracting  bodies  that  the  assumption  of  point  masses 
is  supported.  The  strongest  effect  would  be  that  of  the  Earth  and 
the  largest  term  (J2  term)  to  h<_ve  effect  acts  to  the  fifth 
power  of  the  distances  from  the  Earth  making  this  term  negligible. 
Consider  Mac  Cullagh's  Formula  (Ref  13)  which  gives  the  potential 
for  an  attracting  body  of  any  shape  at  a  distance  from  the  attracting 
body  which  is  large  compared  with  the  body's  over-all  dimensions. 

V  =  f"2  "  Ir3  (A  +  B  +  C  -  31) 

where  I  =  r2  sin  20  dm 
and  A  +  B  +  C  =  2  r2dm 

where  A,  Bf  and  C  are  the  principle  moments  of  inertia  and  I  the 
moment  of  inertia  of  the  body  about  a  line  drawn  from  the  center 
of  gravity  of  the  body  to  the  satellite.  From  this  equation  it 
can  be  seen  that  the  satellite  does  not  have  to  be  very  many 
times  the  radius  of  the  body  away  before  the  second  term  on  the 
right-hand  side  becomes  negligible. 

The  second  assumption  is  obvious  when  you  consider  the 
tremendous  resources  needed  to  orbit  a  satellite  of  any  great 
size.  Even  a  large  space  colony  would  only  have  to  be  ten  miles 
across  to  meet  the  needs  of  the  people  and  one  this  size  would 
exert  no  influence  upon  the  attracting  bodies. 

The  effect  of  the  other  planets  in  the  solar  system  on  the 
Sun  will  be  taken  into  account  by  the  eccentricity  and  formula 
used  for  the  Sun's  motion.  This  will  account  for  most  of  the 
direct  effect  since  direct  effect  upon  the  Earth  and  Moon  will 


not  be  as  great  as  upon  the  Sun.  The  indirect  perturbations  are 
larger  and  tend  to  change  the  origin  to  the  center  of  mass,  or 
barycenter,  of  the  Sun  and  planets,  Schutz  and  Tapley  studied 
the  effect  of  other  planets  using  an  Ephemeris  Model,  Numerical 
integration  was  performed  for  L5  with  a  fixed  step  size  and 
compared  with  results  of  integration  neglecting  the  planets  at 
the  same  step  size.  The  maximum  change  in  position  was  488  km 
which  was  small  enough  to  not  significantly  affect  the  motion 
during  2500  days  considered  but  it  would  have  a  significant  effect 
over  longer  periods.  The  effect  is  slight  enough  that  a  small 
controller  could  offset  its  influence. 

Out-of-plane  considerations  have  been  dealt  with  using 
intrinsic  solutions  before*.  It  is  necessary  to  consider  effects 
due  to  the  inclination  of  the  Moon  and  satellite  with  respect  to  the 
ecliptic.  The  lunar-orbit  plane  is  inclined  to  the  ecliptic  by  5,14 
degrees.  Consequently,  if  the  colony  is  initially  in  a  coplanar 
orbit,  there  are  different  rates  of  regression  of  the  lines  of  nodes 
of  the  Moon  (due  to  solar  perturbations)  and  of  the  satellite  (due 
to  lunisolar  perturbations) .  Although  both  Moon  and  satellite 
maintain  nearly  constant  inclinations  on  the  ecliptic,  their 
orbit  planes  mutually  precess  and,  in  time,  are  mutually  inclined 
by  up  to  10  degrees.  However,  as  Danby  discusses  (Ref  13) ,  the 
angular  momentum  of  the  solar  system  is  almost  totally  coplanar 
giving  the  system  an  invariable  plane.  Any  orbit  within  the 
system  that  leaves  this  plane  will  be  drawn  back  in  and  the 
terms  due  to  perturbation  of  the  Moon's  inclination  are  periodic. 
Heppenheimer  (Ref  20)  estimates  the  optimal  satellite  inclination 
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using  phase  equilibria  and  finds  secular  and  inclination  - 
type  resonances  to  be  small  and  stable. 

Restricted  Four-Body  Equations  of  Motion 

The  general  procedure  followed  to  obtain  real  world  orbits  is 
to  use  current  numerical  integration  schemes  on  planar  restricted 
four-body  equations  of  motion.  In  these  equations,  the  motion  of 
the  Sun  is  taken  as  an  unperturbed  ellipse  with  respect  to  the 
Earth  -  Moon  barycenter  and  the  lunar  motion  is  given  by 


The  motion  of  the  satellite  is  given  by 


(4) 


where  the  subscripts  m,  s,  and  c  stand  for  the  Moon,  Sun,  and 
satellite,  respectively.  These  equations  are  derived  in  Appendix 
A  and  were  used  by  Heppenheimer  (Ref  20).  They  use  a  rectangular, 
nonrotating,  Earth  -  centered  coordinate  system.  These  equations 
reduce  to  the  restricted  three-body  problem  is  the  perturbation  of 
the  fourth  body  (Sun)  is  removed.  The  remaining  terms  are 

r  2  =  (x  -  x  ) 2  +  (y  -  y  ) 2  (Sun  -  Moon  distance) 
ms  ms  ms 
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+  (3)  e  ^  sin  2  n  t.  +  (13)  a  ^  sin  3  n  t  (6) 

- Q  <'  — ,  C*  c- 

4  12 

so  that  the  Sun's  coordinates  are  given  by 

x  =  y  x  +  rp  cos  f 
sms  s 

y  =  U  y  +  rp  sinf 
s  m  s  s 

rp  ~  a  (1  -  e  2) 
s  s  s 

1  +  e  cosf 
s  s 

Lunar  motion  is  initiated  also  on  the  positive  x-axis.  Refer 
to  Figure  1  for  a  diagram  of  the  initial  condition  configuration. 

These  equations  of  motion  arc  cast  into  first-order  form. 

The  initial  conditions  of  the  periodic  stable  orbits  found  by 
the  various  authors  are  to  be  determined  and  transformed  into 
the  preceding  four-body  coordinate  system.  The  equations  of  tire 
four-body  model  are  integrated  using  a  fourth-order  Adams  - 
Bashforth  predictor  and  Adams  -  Moulton  corrector  with  a  fourth- 
order  Ru.nge  -  Kutta  integrator  as  a  starter.  This  entire 
integration  package,  called  ODE,  was  developed  by  Shanpino  and 
Gordan  at  Sandia  Laboratory.  It  automatically  adjusts  the  order 
and  step  size  to  control  the  local  error  per  unit  step  in  a 
generalized  sense.  It  is  the  integration  package  used  throughout 
this  report  on  the  CDC  6600  and  CYBER  tie-in  computers  at  the  Air 
Force  Institute  of  Technology. 

Numerical  experiments  will  not  have  to  be  performed  to 
determine  a  proper  step  size  for  integration  of  the  equations 
over  a  given  range.  ODE  will  optimize  the  step  size  preventing 
large  round-off  error  due  to  a  small  step  size  over  a  large 
number  of  integration  steps  or  truncation  error  due  to  a  large 


20 


1 


stop  size  over  many  steps- 

Other  constants  used  by  Hop;  onhnimt.r  which  will  ho  compared  to 
the  other  authors'  constants 

Unit  mass  -  1  (Earth  +  Moon)  (7) 

Lunar  mass  -  p  =•-  0.01215 

Unit  distance  --  3.8441  >:  10^  m  (Ecirth  -  Moon  distance) 

Unit  time  -  4.3484167  days 

T.U. 

Unit  velocity  -  1023.17  m/scc 

Unit  acceleration  =  0.00273  ra/sec2 
The  unit  time  corresponds  to  one  radian  per  time  unit  mean 
motion  of  the  Moon  for  a  27.321661  day  sidereal  period  ai.d 
0.07480133  radian  per  time  unit  mean  motion  of  the  Sun  for  a 
365.256305  day  sidereal  year. 
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III.  WHEELER'S  ORBIT 
Overview  of  Wheeler's  Work 

In  Wheeler's  study  (Ref  49)  the  equations  of  motion  for  a 
satellite  near  L4  in  a  planar  restricted  four-body  problem  are 
derived  and  tested.  A  computer  algorithm  for  finding  periodic 
motion  is  formulated  and  the  initial  starting  point  of  the  Sun, 

Moon  and  Earth  in  line  is  used.  Rabe's  periodic  orbits  in  the 
three-body  problem  are  used  as  starting  conditions  to  begin 
searching  for  four-body  motion.  He  employs  a  linearization  of 
small  displacements  about  the  non-linear  periodic  orbit.  The 
periodic  orbit  found  has  a  period  in  resonance  with  the  lunar 
synodic  month,  29.5305882  days.  Wheeler  again  calculates  this 
same  orbit  by  slowly  increasing  the  mass  of  the  Sun  from  zero 
to  its  real  value  and  presents  a  synodic  period  orbit  about  L4. 

He  proves  this  orbit  is  generated  mathematically  from  the  L4  point 
and  not  from  Rabe's  orbits.  A  successful  stability  analysis  is 
performed  on  the  orbit  with  the  orbit  found  to  be  stable. 
Assumptions  and  Coordinate  System 

The  assumptions  made  using  wheeler's  very  restricted  four- body 
model  (Ref  49)  include: 

1)  The  Sun,  Earth  and  Moon  are  considered  to  be  point  masses, 

2)  The  mass  of  the  satellite  is  negligible  when  compared  to 
the  other  three  bodies  and  it,  therefore  does  not  affect  their 
motion. 

3)  The  gravitational  effects  of  the  other  planets  in  the 
solar  system  can  be  ignored. 
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4)  The;  motion  of  till  four  bodies  in  limited  to  one-  piano. 

5)  The  Earth  and  Moon  move  in  circular  orbit;;  about  their 
barycenter  at  a  constant,  rotation  rate. 

G)  The  Earth  -  Moon  barycenter  moves  in  a  circular  orbi t 
about,  the  Sun  at.  a  constant  rotation  rate. 

The  difference  of  Wheeler's  model  is  obviously  the  assumption 
of  circular  orbits  instead  of  the  actual  perturbed  orbits  for  the 
Earth,  Moon  and  Sun.  Later  in  tin's  report.  Wheel  or 'f.  orb' t  wil  1. 
be  reproduced  using  the  truth  model  by  removing  the  perturbative 
forces  from  t.hc  Moon  and  Earth  orbit.  These  forces  need  to  be 
accounted  for  to  correctly  model  orbital  motion  in  cislvnar  space. 

Refer  to  Figure  2  for  Wheeler's  rotating  coordinate  system. 
The  satellite  orbiting  about  L4  remains  in  a  rotating  coordinate 
system  with  its  center  at  the  Earth  -  Moon  barycenter  and  the  Moon 
on  the  negative  x-axis.  So  the  coordinate  system  rotates  with  the 
period  of  the  Moon's  synodic  period  and  the  entire  system  rotates 
in  turn  about  the  barycenter  between  the  Earth  -  Moon  barycenter 
and  the  Sun.  Its  rotational,  rate  is  that  of  one  sidereal  year. 
Constants 

Pertinent  constants  used  in  Wheeler's  paper  are 

Solar  mass  =  m  =  328900.12  (8) 

s 

Mean  solar  semimajor  axis  -  a  =  388.82028 

s 

Solar  eccentricity  =  e  =0.0 

s 

Solar  eccentricity  is  equivalent  to  the  movement  of  the 
Earth  -  Moon  barycenter  about  the  Sun. 


Pig  2.  Wheeler's  Four-Body  Configuration 


.0803489351 


Solar  mean  motion  =  n  = 

s 

From  Figure  2,  0  is  normalized  to  1,0  and  so  the  entire  rotation 
of  the  Earth  -  Moon  system  with  respect  to  an  inertial  frame 
w  =  0  +  a 

w  =  1.0  +  .0808489351  =  1.0808489351  (9) 

which  is  equivalent  to  lunar  mean  motion. 

If  the  mass  of  the  Moon  is  set  to  some  arbitrary  valve,  p, 

and  the  sum  of  the  masses  of  the  Earth  and  Moon  is  set  equal  to 

► 

1,  then  the  masses  of  the  Earth  equals  1  -  p.  These  values  also 
correspond  to  the  normalized  Earth  -  Moon  distance  between  their 
barycenter.  In  Wheeler's  study 

p  =  .0121396054 
Conversion  to  Truth  Model  * 

The  constants  need  to  be  converted  to  the  truth  model  situation 
and  the  initial  conditions  of  Wheeler's  orbit  which  will  be  used 
to  generate  the  truth  model  orbit  must  be  transformed  into  the 
coordinate  system  used  for  the  truth  model. 

The  initial  conditions  of  Wheeler's  1/1  resonance  orbit  are 
(Ref  49)  x  =  -  0.72418782459  x  =  0.07948061949  (10) 

y  =  0 . 8x^68639689  y  =  0.22438007788 
They  were  generated  from  the  initial  conditions  at  the  L4  point 
x  =  -  0.4878603946  x  =  0.0  (11) 

y  =  0.8660254038  y  =  0.0 

for  n»s  =  0.0  (three-body  problem)  and  increasing  up  to  its  present 
value  in  a  very  restricted  four-body  problem.  The  initial  con- 
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ditions  (10)  provide  periodic  motion  in  the  restricted  four-body 


model  for  periods  which  are  integer  multiples  of  the  lunar  synodic- 
month,  Figure  3  describes  Wheeler's  initial  condition  configuration 
with  the  Sun  and  Moon  starting  on  the  negative  x-axis  of  the 
rotating  B  coordinate  svstem  with  angles  0  and  a  equal  zero. 

Wheeler's  initial  conditions  must  bo  transformed  from  a 
rotating  Earth  -  Moon  barycentcr  frame  to  an  Earth  -  centered 
inertial  frame.  Referring  to  Figure  4,  the  frames  are  initially 
lined  up  and  the  distance  from  the  origin  of  the  truth  mode]  E 
to  the  initial  starting  point  is 

-P/E  =  — P/B  +  — B/E 

where 


—P/B  =  x  £i  +  ^2  and  —B/E  =  -  1J  % 


combining 


-P/E  =  (X  -  P)  %  +  y  —2 


(12) 


in  terms  of  Wheeler's  coordinates.  For  the  velocity  vector 
e  e  e 

v  =  V  +  V 
-P/E  -P/D  -B/E 

e  ea  e  c  ca  e 

=  v  ,  +  w  xr,  +  v  .  +  w  x  r  , 

-P/B  -  -P/B  -B/E  -  -B/E 

where 


— P/B  =  X  +  y  ^2 


—B/E  =  ^  — 1  =  ° 


ea  ' 
w  =  0  e 
~  -3 
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e, 
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4,  Transformation  From  Wheeler  Frame 


To  Truth  Frame 


substituting 


e 

v 

-P/E 


v  • 

=  0  e3  x  (-  p  e^)  +  x  e^  +  y  e_2  +  0  e 3  x  (x  e  +  y  e^) 


Zp/E  =  (x  -  0  y)  £3  +  (y  +  0  x  -  0  p)  e2 

Equations  (12)  and  (13)  are  in  terms  of  Wheeler's  coordinates 
and  must  be  rotated  into  truth  model  coordinates.  The  rotation 
matrix  is 

a  .  ae  e 
^P/E  ~  — P/E 


— P/E 


cos  0  -  sin  0 
sin  0  cos  0 


—P/E 


Multiply  the  rotation  matrix  times  equations  (12)  and  (13)  to 
get  the  coordinates  in  terms  of  the  truth  model  frame. 


x  =  (x  -  p)  ftos  0  -  y  sin  0 

cl 


(14) 


y  =  (x  -  p)  sin  0  +  y  cos  0 

a 

•  •  • 
x  =  (x  -  0y)  cos  0  -  (y  +  0  x  -  0  p)  sin  0 

cl 

y  =  (x  -  0y)  sin  0  +  (y  +  0  x  -  0  p)  cos  0 

cl 

Since  Heppenheimer ' s  model  initially  positions  the  Sun  and  Moon 
on  the  positive  x-axis,  the  Sun  and  Moon  starting  positions  will 
be  rotated  180°  to  the  negative  x-axis  and  0=0°.  In  Wheeler's 
frame  0=1  but  there  is  a  time  difference  between  the  two  frames 
which  will  affect  the  transformation  of  the  velocity  initial 
conditions.  In  Wheeler's  frame  the  period  of  rotation  is  one 
lunar  synodic  month,  29.5305882  days,  whereas  in  the  truth 
model  frame  the  period  of  rotation  is  one  lunar  sidereal  month, 
27.32166101  days.  So  the  slower  truth  model  frame  must  have  the 
Wheeler  velocities  scaled  down  by  a  factor 
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27.32166101  _  1 

29.5305882  =  1,0808489351 

So  the  transformation  equations  (14)  are  now 

x  =  x  -  y  (15) 

d 

ya  =  y 

x  =  -  y  +  x/1. 0808489351 

cL 

y  =  x  -  y  +  y/1. 0808489351 

d 

Substituting  Wheeler's  initial  conditions  (10)  we  obtained  initial 
conditions  for  the  truth  model 

x  =  -  0.73632742999  x  =-0.74215104336  (15) 

y  =  0.81568639689  y  ="0.52873127999 

Wheeler's  constant  (8)  will  be  used  in  the  truth  model  along 

t 

with  his  value  of  y  and  the  transformed  initial  conditions.  The 
solar  mean  motion  of  Wheeler's  model  must  be  converted  to  the 
slower  truth  model  frame. 


%  =  0.0808489351 
w  1.0808489351 


0.7480133 


(17) 


The  true  anomoly  of  the  equation  of  the  center  (6)  will  be 


rotated  to  the  negative  x-axis  by  having  H  or  180°  added  to  it. 
The  Moon  initial  conditions  will  also  begin  on  the  negative 
x-axis  to  correspond  with  Wheeler's  scenario. 


Verification  of  the  Truth  Model 

The  algorithm  for  integration  of  the  equations  of  motion 
and  plotting  the  orbits  formulated  is  given  in  Appendix  C. 
Wheeler's  assumptions  are  applied  to  the  truth  model  to  verify 
that  the  same  orbit  is  reproduced.  Wheeler's  orbit  of  1/1 
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resonance  is  shown  in  Figure  5.  Since  Wheeler's  assumptions  involve 

circular  orbits  the  eccentricity  of  the  Sun  is  zero  ( e -  0.0) 

and  the  Moon's  initial  conditions  are 

x  =  1.0  x  =  0.0  (18) 

m  m 


y 


m 


0.0 


y  =  Hr  =  -  1,0  where  u  =  (Gm  +  m  ) 
m  yp  e  m 

=  1.0 


The  truth  model  equations  of  motion  for  lunar  motion  are  simplified 
by  eliminating  the  perturbative  effect  of  the  Sun.  Equations 
(1)  and  (2)  now  become 


x  +  x  =0 
m  m 


m 


y  +  y  =0 


(19) 


m 

The  equations  for  satellite  motion  (3)  and  (4)  remain  the  same. 

A  period  of  rotation  of  the  truth  model  is  2V  or  6.2831853 
radians  or  time  units.  But  Wheeler's  period  occurs  when  the 
bodies  (Sun,  Earth,  Moon)  are  all  lined  up  again  which  is  a 
longer  period  of  time 

__  29.5305882  days  r  .  .  . . 

211  x  ~  ■ 1  =  6.791174148  time  units  (20) 

27.32166101  days 


so  the  orbit  period  will  be  considered  this  number  of  time  units 
and  the  orbit  found  will  be  transformed  in  wheeler's  frame  for 
comparison  and  plotting  purposes. 

The  equations  of  motion  were  integrated  for  over  fifty 
years  and  the  orbit  found  (Figure  5)  is  identical  to  Wheeler's. 
Figure  6  and  7  show  three  and  20  orbits  respectively.  Points 
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Fig  5*  Wheeler's  Very  Restricted  Four-Body 

1/1  Periodic  Orbit,  1  Orbit 
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Fip  6.  Wheeler's  Very  Restricted  Four-Body 

1/1  Periodic  Crbit,  3  Crbits 
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along  the  periodic  stable  orbit  were  matched  identically  over 
the  same  time  period  less  computer  round-off  in  the  fifth  and 
smaller  digits.  The  orbit  is  in  the  same  orientation  about 
L4  with  an  approximate  semimajor  axis  of  80,000  miles  and 
semiminor  axis  of  40,000  miles.  This  is  strong  evidence  that 
Wheeler's  equations  for  the  very  restricted  four-body  model  are 
correct  and  that  the  truth  model  from  Heppenheimer  is  correct 
and  can  be  reduced  to  the  three-body  problem.  The  next  step  is 
to  apply  the  full  model  to  wheeler's  orbit. 

Truth  Model  Application 

It  greatly  simplifies  the  four-body  problem  to  start  out 
with  the  Sun,  Moon  and  Earth  lying  on  a  straight  line.  Thus, 
in  order  for  periodic  motion* to  occur,  the  equations  must  be 
integrated  forward  to  the  exact  time  when  these  initial  con¬ 
ditions  occur  again  and  when  the  satellite  has,  of  course, 
returned  again  to  its  original  position  about  L4.  This  orbit 
being  examined  returns  after  one  lunar  synodic  month  in  the 
very  restricted  four-body  case.  But  as  the  perturbations  increase 
by  adding  more  segments  to  the  truth  model  the  orbit  will  no 
longer  return  to  exactly  the  same  spot  or  at  exactly  the  proper 
time.  The  synodic  period  will  continue  to  be  used  as  a  basis 
for  looking  at  orbit  positions  to  see  how  they  are  changed  by 
the  more  realistic  forces  of  the  four-body  problem. 

Perturbations  were  added  one  at  a  time  to  observe  the 
effects  upon  stability  and  periodicity  noting  that  one  effect 
could  offset  part  of  another  one  when  all  the  forces  are  inter- 
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Fi*  8.  Wheeler  Orbit  With  Solar  Eccentricity, 


1  Orbit 
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Fig  9*  Wheeler  Orbit  With  Solar  Eccentricity , 


3  Orbits 
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0.72 


Fig  10.  Wheeler  Orbit  With  Solar  Eccentricity, 


20  Orbits 
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Fig  11.  Wheeler  Orbit  With  Sun-Perturbed  Lunar 


Motion,  1  Orbit 
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Wheeler  Orbit  With  Sun-Perturbed  Lunar 

Motion,  3  Orbits 
41 
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acting.  The  eccentricity  of  the  Sun  was  input  and  the  orbit 
integrated  for  ten  years.  Figures  8,  9  and  10  depict  the  orbit 
for  one,  three  and  twenty  periods  respectively.  Note  there  is 
only  a  slight  perturbation  away  at  each  position  compared  to  the 
previous  orbit.  The  change-  from  one  orbit  to  the  next  is  at 
most  650  km  anywhere  in  the  orbit  and  the  basic  orbit  shape  or 
position  about  L4  does  not  change  after  ten  years.  Orbits  after 
ten  years  are  no  more  than  1000  km  away  from  the  first  orbit  and 
tend  to  return  to  the  starting  orbit. 

Next  the  perturbative  effect  of  the  Sun  upon  lunar  motion 
is  returned  to  the  problem  by  using  equations  (1)  and  (2) .  The 
orbit  is  now  highly  perturbed  changing  shape  somewhat  and  up  to 
4000  km  off  initial  orbit  position  the  second  time  around.  The 
orbit  slowly  begins  to  drift  out  of  the  quadrant  of  Wheeler's  frame 
after  only  thirteen  months.  This  is  not  so  damaging,  but  the 
changing  shape  of  the  orbit  shows  little  consistency  in  the 
orbit,  as  Figures  11,  12  and,  especially,  13  will  attest.  After 
ten  months  the  orbit  approaches  the  Moon  close  enough  to  be 
thrown  out  of  the  Earth  -  Moon  system  completely. 

The  truth  model  appeared  to  prove  the  Wheeler  orbit 
unstable.  Analysis  was  made  of  the  data  to  try  to  determine 
the  cause  of  the  instability.  The  eccentricity  of  the  Sun  is 
removed  with  little  effect  on  the  orbit,  showing  that  the 
direct  influence  of  the  perturbing  Sun  does  not  change  the 
stability  of  a  periodic  orbit  with  its  small  deviations.  The 
two  terms  of  the  Sun's  effect  on  lunar  motion  were  analyzed 
and  found  to  offset  each  other  to  such  a  degree  that  removing 
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one  or  the  other  ear.  s  the  orbit  to  deteriorate  even  more 
rapidly.  Finally,  the  phase  difference  between  the  Sun  and 
Moon  is  examined  in  the  truth  model.  This  is  thought  to  be 
the  primary  cause  for  the  loss  of  stability.  A  program  is 
written  to  minimize  the  difference  in  the  phase  of  the  two 
bodies  by  adjusting  the  solar  mean  motion.  If  the  Sun  and 
Moon  were  held  in  phase,  this  author  believes  the  stability 
would  return.  The  optimization  was  not  very  efficient  but 
after  three  cycles  had  improved  stability  50t.  The  orbit  is 
indeed  sensitive  to  the  Moon  getting  out  of  phase  with  the  Sun. 

During  this  phase  analysis  an  oversight  came  to  light. 
Although  the  lunar  motion  now  includes  perturbing  forces  of 
the  Sun  to  force  the  Moon  out  of  circular  motion,  the  lunar 
initial  conditions  arc  still  the  circular  orbit  conditions 
determined  earlier  (IS) .  Elliptical  initial  conditions  must 
be  found.  Kolcnkiewi.cz  and  Carpenter  (Ref  26)  developed  a 
model  of  the  Moon's  orbit  using  a  three-body  system  which  is 
explained  in  detail  in  the  Chapter  Four.  Since  a  perturbed  orbit 
was  obtained  using  the  two  primary  bodies.  Earth  and  Sun,  which 
affect  it,  this  orbit  is  considered  accurate  enough  to  supply 
initial  conditions. 

x  =  -  0 . 99220573479  x  =  0.0  (21) 

m  m 

y  =  0.0  y  =  1.00990709043 
m  m 

When  these  perturbed  lunar  initial  conditions  were  first  input, 
calculator  accuracy  to  eight  places  were  used. 

x  =  0.99220574  x  =  0.0  (22) 

m  m 
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y  =  0.0  y  =  -1.01055416 
m  m 

Now  the  Moon  and  Sun  are  more  in  phase  and  the  orbit  is  reasonably 
stable.  The  orbit  varies  by  no  more  than  5%  over  the  first  24 
months  and  the  orbit  shape  or  size  does  not  change.  It  varies 
by  no  more  than  2%  from  orbit  to  orbit.  Refer  the  Figures  14, 

15  and  16,  for  plots  of  the  orbit.  The  equations  of  motion  are 
integrated  for  over  100  years  and  the  orbit  moves  slowly  around 
the  Earth  in  a  clockwise  direction  away  from  the  Moon,  returning 
to  the  original  orbit  about  L4  after  94  years  with  agreement  to 
within  three  decimal  places.  Refer  to  Figure  20  which  plots  one 
point  on  every  fourth  orbit  for  100  years.  It  is  sixteen  years 
before  the  orbit  begins  to  leave  this  quadrant. 

Even  more  surprising  re'sults  occured  when  the  lunar  initial 
conditions  from  the  computer  (21)  are  input.  Although  these 
initial  conditions  differ  only  in  the  fourth  decimal  place,  the 
orbit  (Figure  17,  18  and  19)  now  moves  slowly  in  a  counter¬ 
clockwise  direction  toward  the  Moon  and  would  again  return  to  the 
same  orbit  about  L4  after  57  years.  The  orbit  begins  to  leave 
the  quadrant  after  ten  years  but  both  of  these  orbits  will  move 
into  the  Moon's  vicinity  and  be  thrown  out  of  cislunar  space. 

The  apparent  drift  of  the  orbit  is  extremely  consistent.  The 
Wheeler  frame  for  plotting  the  orbits  is  investigated  to  dis¬ 
cover  some  reason  for  the  apparent  drift.  The  Wheeler  frame 
rotates  at  a  constant  rotation  rate  to  coincide  with  the  lunar 
synodic  month.  The  elliptical  motion  of  the  Moon  and  the 
variations  in  its  rotation  rate  is  causing  the  reference  frame 
to  drift  rather  than  the  orbit.  A  new  algorithm  is  formulated 
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Fig  14.  Truth  Model  of  Wheeler  Orbit,  Calculator 

Accurate  Lunar  Init.  Cond.,  1  Orbit 
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Fig  15-  Truth  Model  of  Wheeler  Orbit,  Calculator 

Accurate  Lunar  Init.  Cond.,  3  Orbits 
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Fig  16.  Truth  Model  Of  Wheeler  Orbit,  Calculator 

Accurate  Lunar  Init.  Cond 2®  Orbits 
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Fig  17.  Truth  Model  of  Wheeler  Orbit,  Computer 

Accurate  Lunar  Init.  Cond.,  1  Orbit 
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Fig  18.  Truth  Model  of  Wheeler  Orbit,  Computer 

Accurate  Lunar  Init.  Cond.,  3  Orbits 
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Fig  20.  Drift  of  Wheeler  Orbit  Over  60  Years 
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to  lock  the  reference  frame  to  the-  Moon  wherever  it  is  and  keep 
it  on  the  negative  x-axis.  When  this  is  done  the  Wheeler  truth 
model  is  found  to  be  very  stable  with  less  than  a  5%  change  in 
position  over  45  months  and  10%  over  seven  years.  It  changes 
from  orbit  to  orbit  no  more  than  4000  km.  The  orbit  slowly 
wanders  back  and  forth  for  over  100  years  and  never  leaves  the 
L4  quadrant.  The  orbit  slowly  changes  shape  over  the  years 
depending  on  its  position  in  the  quadrant  in  relation  to  the 
other  bodies. 
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IV.  KOLENKIEWICZ  AND  CARPENTER'S  ORBIT 


Overview 

Kolenkiewicz  and  Carpenter  (hereafter  termed  K  &  C)  use  a 
restricted  four-body  model  in  which  the  three  principle  bodies  are 
periodic,  coplanar,  and  obey  the  equations  of  motion  and  have  no 
mean  orbital  eccentricities  (Ref  26) .  They  investigated  by 
trigonometric  series  the  possibility  of  a  coplanar  monthly 
periodic  motion  in  the  general  vicinity  of  L4.  This  investigation 
of  the  Sun  perturbed  Earth  -  Moon  triangular  point  yields,  in 
addition  to  a  small  unstable  orbit,  two  similar  but  not  identical 
stable  periodic  orbits  about  50%  larger  than  Schechter's  stable 
orbits,  one  synchronized  with  the  Sun  in  its  motion  around  L4; 
the  other  180°  out-of-phase!  The  two  stable  orbits  are  periodic 
with  respect  to  the  synodic  system,  make  one  loop  about  the 
triangular  point,  and  are  elliptical  in  shape  (see  Figure  21) . 

The  orbits  have  an  approximate  semimajor  axis  of  90,000  miles 
and  a  semiminor  axis  of  44,000  miles.  The  major  axis  is 
perpendicular  to  the  line  joining  the  Earth  and  L4.  The  particle 
describing  the  orbits  is  synchronized  with  the  Sun  so  their 
angular  positions  almost  coincide  when  the  particle  crosses  one 
of  the  axes  of  the  ellipse. 

Assumptions  and  Coordinate  System 

The  assumptions  made  using  K  &  C's  restricted  four- body 
model  include: 

1)  Sun,  Earth  and  Moon  are  considered  point  masses, 

2)  The  gravitational  effects  of  other  planets  are  ignored. 

3)  Only  gravitational  forces  are  considered. 
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Fig  21.  Kolenkiewicz  and  Carpenter's  Results 
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4)  The  motion  of  all  four  bodies  is  limited  to  one  plane. 

5)  The  Moon's  orbit  is  perturbed  by  the  Sun  and  Earth  alone. 

6)  The  Sun's  orbit  is  perturbed  by  the  Moon  and  Earth  alone. 

The  only  difference  from  the  truth  model  is  the  unrealistic 

orbit  obtained  for  the  Sun  which  gives  it  an  approximate 

eccentricity  of  only  .0007313  instead  of  the  actual  value  of 

.0168.  This  should  have  little  effect  since  removing  the  solar 

eccentricity  from  the  truth  model  on  Wheeler's  orbit  gives  a 

slightly  larger  orbit  and  a  bit  legs  stability. 

The  coordinate  system  (Figure  21)  is  the  same  as  the  truth 

model,  an  Earth  -  centered  inertial  system  with  the  Moon  and  Sun 

lying  along  the  positive  x-axis  initially  and  the  triangular 

point  lying  60°  off  that  akis.  There  is  another  internal  coordinate 

system  for  locating  the  satellite  in  relation  to  the  triangular 

point  rotating  in  advance  of  the  Moon,  Its  origin  is  located  at 

<t>  =  60°;  r  =  a  =  3.831841237  x  108m  (23) 

— o 

from 

n,2  a3  =  Gm 
3  e 

Constants 

Many  of  the  constants  are  in  the  SI  system  and  will  remain 

or  be  converted  into  SI  until  calculated  by  the  computer  but 

values  corresponding  to  units  previously  used  will  appear  in 

parenthesis  where  applicable. 

mass  ratios:  m  =  . 

m  1_ 

m  81.30 
e 

m 

—  =  332958.087932061 
m 
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therefore 

Solar  mass  -  in  =  328912.42  (24) 

s 

Vi  =  .01215067 

also 

[ 

Mean  solar  scraimajor  £txis  ~  a  =  149600  >:  10^  ra  \ 

I 

i 

(=  388.8235  mean  E  -  M  distances)  j 

Mean  lunar  scinimajor  ?ixis  -  a  ~  3.847487965  x  10^'  m 

m 

{  =  1.0  mean  E  -  M  distance) 

Solar  moan  motion  =  n  =  129597742" . 38 

3 - — - - — ■ - 

Julian  Century 

=  2,0  x  10  7  radians 
sec 

(~  .08084893  radians  ) 
time  unit 

Solar  moan  motion  =  n  =  173 2559353 ■ " 56 

Julian  Ccnt.ury 

=  2.66  x  10  6  radians 
sec 

(=  1.0808489  radians  ) 
time  unit 

The  mean  motions  correspond  to  -1)10  lunar  synodic  month. 

Geocentric  gravitational  constant  = 

Gm  =  398603  x  10^  mVsec2 
e 

since 

m  =  .98784933 
e 

then 

G  =  4.0350283  x  1014 
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Conversion  to  Truth  Model 

K  &  C's  procedure  used  in  obtaining  a  solution  for  the 
Earth  -  Moon  -  Sun  model  as  well  as  for  the  periodic  triangular 
point  orbits  is  based  on  Musen’s  (Ref  30)  method  with  the  per¬ 
turbations  represented  in  trigonometric  series  with  numerical 
coefficients.  The  solutions  of  the  equations  are  given  in  the 
following  form  (see  Figure  21) 

r  =  (1  =  a)  r  +  8  w  (25) 

—  — o  — 

where  a  and  8  are  the  components  of  the  perturbations,  r^  is  the 
position  vector  in  a  fixed  reference  ellipse,  and 

<i)  =  (1/n)  (dr  /dt) 


The  mean  motion,  n,  is  given  but  can  be  calculated  as  a  check 
from  Kepler's  law  which  is  n2  a2  =  u2  where  a  is  the  semimajor 
axis  of  the  reference  ellipse.  All  the  reference  ellipses  used 
have  zero  eccentricity  (circles)  so  rQ  =  a. 

The  functions  a  and  8  are  represented  by  the  trigonometric 
series 


=  £ 
k=0 


(a,  coskS  +  a  sin  k9) 


(26) 


where 


0 


£ 

k=0 


(8^°  cosk9  +  8^  sin  k8) 


(27) 


0  =  (n  -  n  )  t. 
m  s 

The  basic  inertial  coordinate  system  matches  the  truth  model 
and  needs  no  transformation.  However,  the  initial  conditions  need 
to  be  determined  from  the  orbits  of  the  Sun,  Moon  and  satellite 


and  these  are  described  using  constant  coefficients  of  a 


trigonometric  series  (Tables  1  -  4) ,  The  initial  conditions  of  each 
orbit  can  be  determined  from  the  coefficients  at  t  =  o. 


Expanding  on  equation  (25)  and  stating  the  velocity  equation,  r 


r  =  (1  +  a)  r  + 
—  — o 


n  V  r 


r  =  (<*  -  8.  M2  ]  r^  +  I&  +  1  +  a]  AT2  a„ 

n  V  r 


(28) 

(29) 


n  r 


o 

where  r  is  in  the  a,  direction  in  the  truth  model  frame.  The 
—o  —1 

values  of  a  and  8  can  be  found  from  equations  (26)  and  (27)  for 
t  =  o  and  9=0°  for  the  Sun  and  Moon  and  0=60°  for  the  satellite. 
The  equations  for  a  and  8  are 


a  =  (n 

m 


8  = 


(n 


m 


n  )  E  (-  k  a  c  sin  k0  +  k  a  S  cosk9) 

5  .  «  iC  <K 

k=0 


n  )  E  (-  k  8.  C  sin  k0  +  k  a  S  cosk0) 

S  K  JC 

k=0 


(30) 

(31) 


Solutions  for  the  motion  of  the  Moon  and  Sun  are  found  using 
K  &  C's  same  four-body  equations.  Starting  with  the  Sun  constrained 
to  move  in  a  circular,  coplanar,  Keplerian  orbit  with  respect  to 
the  Earth,  the  equations  of  motion  are  solved.  Trigonometric 
coefficients  a  and  8  describing  the  Moon's  perturbed  orbit  are 
thus  obtained.  The  role  of  the  bodies  is  reversed,  the  Moon's 
motion  is  constrained  to  move  in  the  perturbed  orbit  defined  by 
a  and  8,  and  the  equations  of  motion  for  the  Sun  are  solved.  The 
a  and  8,  coefficients  describing  the  Sun's  perturbed  orbit  are  thus 
obtained.  The  roles  of  the  bodies  are  reversed  again  and  again, 
each  time  using  the  latest  acquired  a  and  8  coefficients  of  each. 
Ultimately  the  values  of  a  and  8  converge  for  each  body.  The 
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Table  1  Three-body  sun  solution,  <?  «  0* 


k 

«*<•>  X  10* 

3t«"  X  10* 

0 

0  045105 

0.000000 

1 

30.949S6S 

31.492851 

2 

-0.004947 

-0.005012 

3 

0.047329 

0.047319 

4 

-0.000005 

-0.000005 

5 

0.000184 

0.000183 

6 

0.000000 

0.000000 

7 

0.000001 

0.000001 

Table  2  Three-bod  y  moan  solution. 


«»<*'  X  10* 


3."'  X  10* 


0 

1 

2 

3 

4 

5 

6 
t 
8 
9 

10 


-906.915740 
287.606767 
— 7173.506S63 
-7.507078 
6.02S443 
-0  003392 
0.032454 
0  000011 
0.000187 
0.000000 
0.000001 


0  000000 
-609  076345 
10202.254541 
7.212259 
5.719334 
0  005S16 
0.027566 
0  000025 
0.000163 
0  000000 
0.000001 


Table  3  Periodic  orbit  I,  *  ■=  60* 


k 

«*<•'  X  10* 

a*‘*>  X  10* 

d»<*>  x  io* 

/»*«•»  X  10* 

0 

- 19171. 00S123 

0.000000 

74753.542768 

0.000000 

1 

1S7S01 . 13597S 

1717S. 314016 

-13120  769748 

-377986  16.5215 

•> 

11131  030003 

—  3722  87205s 

2352.545921  ' 

lS027.0lo46-l 

3 

—2874.47241* 

737.o(fc02S 

—  637 . 564775 

-2521 . 769547 

4 

582. 134781 

— 176.9SS257 

173 .574850 

518.179570 

5 

-123.327707 

47.337357 

-46.872029 

-110.607813 

6 

26  370S43 

-12.692022 

12.244230 

24.010789 

7 

— 5.8300S3 

3.27140S 

-3.092692 

—5.355760 

S 

1.327555 

-0.827212 

0.781757 

1.245S72 

9 

-0.314593 

0.214693 

-0.205945 

-0.299764 

10 

0.075554 

— 0.05S253 

0.056G69 

0.072270 

11 

— 0.017S29 

0.016064 

— 0.015GS0 

-0.017002 

12 

0.004114 

-0.004335 

0.004214 

0.003918 

13 

—  0.000955 

0  001136 

-0. 001100 

-0.000914 

14 

0.000230 

-0.000296 

0.0002S8 

0.000222 

15 

-0.000057 

O.OOOOSO 

-0.000078 

-0  000055 

16 

0.000014 

-0.000022 

0.000022 

0.000013 

17 

-0.000003 

0.000006 

-0.000006 

-0.000003 

18 

0.000001 

-0.000002 

0.000002 

0.000001 

Table  4  Periodic  orbit  II,  «  —  60* 


k 

at<*>  X  10* 

<.*<•>  X  10* 

0*“>  X  10* 

X  10* 

0 

-18160.912624 

0.000000 

72212. 6SS9S8 

0  000000 

1 

-183627.357659 

- 16818. 0SS205 

11392.917179 

370250.263893 

2 

10460. 90070S 

-3371.790401 

2116.289582 

17696. 741S94 

3 

2715.813738 

-648.117571 

563.594410 

2418.999900 

4 

344.831763 

-152.630065 

151 . 103055 

487.952665 

5 

114.117818 

-40.353414 

40.270023 

102.700732 

6 

24.309814 

-10.697742 

10.365541 

22.014409 

7 

3.274431 

-2.716444 

2.571492 

4.851307 

8 

1.187793 

—  0.674485 

0.637277 

1.115552 

9 

0.278619 

-0.171971 

0. 164972 

0.265686 

10 

0.066360 

-0.046021 

0.044817 

0.063541 

a 

0.013564 

-0.012554 

0.012271 

0.014861  . 

12 

0.003376 

-0.003350 

0.003259 

0.003409 

13 

0.000826 

— 0.000S65 

0  000833 

0.000792 

14 

0.000198 

-0.000222 

0.000213 

0.000192 

13 

0.000049 

-0.000059 

0.000057 

0.000048 

16 

0.000012 

-0.000016 

0.000016 

0.000012 

17 

0  000003 

-0.000004 

0  000004 

0.000003 

18 

0  000001 

-0.000001 

0.000001 

0.000001 

i 
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initial  conditions  of  the  Sun  formed  here  will  not  be  used  because 

the  computed  approximate  eccentricity  is  . 0007313  which  is  not  close 

to  the  actual  ,0168  (Ref  20)  since  other  forces  obviously  act  on 

the  Sun  to  perturb  its  orbit  besides  the  Earth  and  Moon. 

The  lunar  initial  conditions  (21)  are  computed  using  the 

coefficients  from  Table  2,  equations  (26)  and  (31) ,  r  =  a 

o  m 


and 

p  2  =  n  2  a  ^=4.0299  x  1014  m3/sec2 
m  mm 

Since  p  2  =  G  (m  +  m  )  also,  this  was  used  as  a  check  using  values 
m  me 

from  (24) . 

The  satellite  initial  conditions  for  both  orbits  are  computed 

in  a  similar  manner  using  the  coefficients  from  Tables  3  and  4, 

equations  (26),  (27),  (30)  and  (31),  rQ  =  a  from  (23),  9=60°, 

and  p  2  =  Gm  =G(m  +m)  from  (24)  and  the  direction  of  r  and 
c  e  e  c  -o 

m  are 

r  =  a  (cos  9  a,  +  sin  9  a_) 
o  —1  —2 

a j  =  (1/n  )  /  p2  (-  sin  9  a,  +  cos  9  a_) 

—  m  \  /  —  — 1  — z 

V  ro 

The  initial  conditions  of  Orbit  I  which  is  synchronized  with 
the  Sun  in  its  orbit  about  L4  are 

x  =  .7504475430226  x  =  -  .7370404960431  (32) 

y  =  .8127310078686  y  =  .5259745433478 

which  is  the  same  as  Wheeler's  transformed  conditions  (16)  to  two 
places  except  the  x  value  which  is  about  one-hundreth  off. 

The  initial  conditions  of  Orbit  II  which  is  180*  out  of 


phase  are 


x  =  .07542634277321 


x  =-1,003173415668 


(33) 


y  =  .9507111102549  y  =  .254983014774 
Figure  21  shows  the  plotted  coefficients  of  the  Sun,  Moon,  and 
the  two  orbits  compared  to  Schechter's  Reference  X  orbit  (Ref  36) , 
Truth  Model  Application 

Both  orbits  are  run  with  the  initial  conditions  and  constants 
stated  for  over  100  years.  Orbit  I  (Figures  22,  23  and  24)  begins 
with  an  orbit  the  same  shape  but  about  5%  larger  than  the  Wheeler 
orbit.  The  orbit  once  again  changes  shape  very  slowly  and  to  a 
slight  degree  as  it  moves  around  the  Earth  due  to  its  orientation 
to  the  body.  But  this  orbit  drifts  faster  than  Wheeler's  orbit. 

The  orbit  moves  approximately  35,000  km  per  revolution  (month)  and 
begins  to  leave  to  quadrant  after  17  months.  Figure  25  shows  a 
particular  point  on  the  orbit  each  revolution  for  34  years  if  the 
satellite  is  not  thrown  out  of  the  system  by  a  close  approach  to 
the  Moon.  Once  again  this  apparent  orbit  drift  is  deceiving 
because  part  of  it  is  caused  by  the  Moon  drifting  off  the  x-axis 
of  the  Wheeler  reference  frame.  When  the  problem  is  corrected,  the 
orbit  continues  to  drift  but  much  more  slowly,  averaging  less  than 
20,000  km  per  revolution  and  remaining  in  the  quadrant  for  seven 
years.  Within  a  year  after  leaving  the  quadrant  the  satellite  is 
thrown  out  of  the  system  by  a  close  approach  to  the  Moon. 

Orbit  IX  is  integrated  for  over  100  years  (Figures  26,  27 
and  28).  The  first  revolution  is  the  same  size  as  Wheeler's 
truth  model  orbit  but  it  does  not  close  on  itself,  missing  by 
10,000  km.  The  orbit  begins  to  depart  the  Wheeler  frame  L4 
quadrant  at  about  the  same  time  as  Orbit  I,  but  unlike  the  first 
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Pig  22.  Truth  Model  of  K&C  Orbit  I,  1  Orbit 
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Fig  25.  Drift  of  K&C  Orbit  I  Over  60  Years 
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cP.56 


Fig  26.  Truth  Model  of  K&C  Orbit  II,  1  Orbit 
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28.  Truth  Model  of  K*C  Orbit  II,  20  Orbits 
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orbit,  Orbit  II  rapidly  leaves  cislunar  space  after  that.  The 
same  correction  for  the  wheeler  frame  drift  was  made  and  the  drift 
rate  of  the  orbit  once  again  is  found  to  be  slower  than  anticipated 
The  orbit  drifts  an  average  of  30,000  km  per  revolution  and  now 
doesn't  leave  the  quadrant  for  over  10  years  instead  of  an  apparent 
19  months.  Orbit  II  drifts  faster  than  Orbit  I  but  still  has  the 
property  of  being  relatively  stable.  The  shape  of  the  orbit 
indicates  the  possibility  that  with  slightly  different  initial 
conditions  the  orbit  might  come  much  closer  to  closing  upon  itself 
and,  thus,  more  stable.  It  is  interesting  to  note  that  when  both 
orbits  initial  conditions  are  entered  into  the  VRFB  problem  that 
neither  is  periodic  or  very  stable. 
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V.  BARKHAM,  MODI,  AND  SOUDACK'S  ORBIT 

Overview 

Barkham,  Modi  and  Soudack  (hereafter  B,  M  and  S)  offer  a 
theoretical  solution  to  a  four-body  problem  (Ref  1) .  In  their 
model  a  small  particle  moves  in  the  vicinity  of  two  masses,  forming 
a  close  binary,  in  orbit  about  a  distant  mass.  Unique,  uniformly 
valid  solutions  of  this  four-body  problem  are  found  for  motion  near 
both  equilateral  triangle  points  of  the  binary  system  in  terms 
of  a  small  parameter  5,  where  the  primaries  move  in  accordance 
with  a  uniformly  -  valid  three-body  solution.  Accuracy  is 
maintained  within  a  constant  error  0  (?8)  and  the  solutions  are 
uniformly  valid  as  £  tends  to  zero  for  time  intervals  0(£  3) . 

Orbital  position  errors  near  L4  and  L5  of  the  Earth  -  Moon  system 
are  found  to  be  less  than  5%  when  numerically  -  generated  periodic 
solutions  are  used  as  a  standard  of  comparison.  Once  again  the 
lunar  synodic  month  is  used  and  a  21J  periodic  stable  orbit  is 
found  which  makes  two  cycles  per  month.  It  is  a  small  orbit  with 
a  semimajor  axis  slightly  over  6000  km.  The  theoretical  orbit  was 
checked  with  numerical  solutions  of  four-body  perturbation 
equations  of  motion,  from  which  the  theoretical  solution  is 
derived,  simplified  to  three-body  and  with  numerical  solutions  of 
the  theoretical  three-body  equations  where  three-body  solutions 
of  the  Earth  and  Moon  orbits  are  used  to  simplify  the  problem. 

The  two  orbits  computed  have  very  little  difference  and  they  differ 


from  the  theoretical  solution  by  no  more  than  5%. 


Assumptions  _and  Coordinate  System 

The  theoretical  model  assumed 

1)  Sun,  Earth  and  Moon  are  considered  point  masses. 

2)  The  gravitational  effects  of  other  planets  are  dropped. 

3)  The  motion  of  all  four  bodies  lies  in  one  plane. 

4)  All  perturbations  on  the  satellite  except  those  caused 
by  the  Sun  neglected. 

5)  The  Sun  orbits  in  circular  motion. 

6)  Distance  to  Sun  is  decreased  from  infinity  to  present 
distance. 

As  the  Sun's  distance  increases  to  infinity  the  problem  reduces 
to  the  three- body  problem  and  the  particle  contracts  to  the 
libration  point.  The  chocks  to  the  solution  try  to  show  that  the 
Sun  causes  the  primary  perturbations  and  the  satellite  is  relatively 
insensitive  to  small  variations  in  the  orbits  of  the  Earth  and 
Moon.  This  may  be  true  but  indirect  effects  of  the  Sun  upon  the 
Earth  and  Moon  may  cause  large  enough  variations  in  the  orbit  to 
affect  long  term  stability.  A  circular  Sun  has  only  a  slight 
effect  on  the  orbit. 

The  coordinate  system  is  identical  to  Wheeler's  system  (see 
Figures  2  and  29)  with  a  rotating  system  with  the  lunar  synodic 
month  about  the  Earth  -  Moon  barycenter.  The  Moon  and  Sun  start 
initially  on  the  negative  £-axis  or  x-axis. 

Constants 

For  the  Earth  -  Moon  Sun  system  the  following  constants  are 
given  (Ref  1) 


\ 

\ 
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Solar  mass  -  m  =  329794.6384 

s 

Mean  solar  semimajor  axis  =  a^  =  389,17242 

Solar  eccentricity  =  0,0 

Solar  mean  motion  =  n  =  .0808489 

s 

Lunar  mean  motion  =  n  =  1.0808489 

m 

p  =  .01215032 
m 

Conversion  to  Truth  Model 

The  solutions  to  the  satellite  and  Earth  orbits  are  given  by 
trigonometric  series  and  formulas  are  given  to  convert  Earth  to 
Moon  coordinates.  To  find  the  initial  conditions  substitute  t=o 

into  the  series  solutions,  initial  conditions  for  the  satellite 

« 

are 


(34) 


£  =  -.48313894  £  =  .01133797 

n  =  .87089898  g  =  .0044541 

and  for  the  Earth  are 

£  =  .01205487  £  =  0 

n  =  0  n  =  .00023703 

converting  to  the  initial  conditions  for  the  Moon  (Ref  1 ) 

£  =  -.98008937  £  =  0 

n  =  0  T)  -  -.019272 

The  same  translation  is  used  on  coordinates  (34)  and  (35)  as  for 
Wheeler's  system 

r  =  (£  -  p)  a  +  n  a2 
r  =  (£  -  eg)  +  (f|  +  0  £  -  9p)  a2 


(35) 


And  the  velocity  components  have  to  be  scaled  down  for  the  time 


difference  in  the  periods  of  rotation  used  so  divide  by 
1.0808489351,  The  satellite  initial  conditions  are 

x  -  .49528930  x  =  .80040910  (36) 

y  =  .87089898  y  =  .49110833 

and  the  Moon  initial  conditions  are 

x  =  -.99223969  x  =  0  (37) 

y  =  0  y  =  1.0100701 

Truth  Model  App li cation 

The  satellite  is  unstable  when  the  truth  model  is  applied  to 
it.  It  slowly  begins  to  orbit  L4  and  after  half  a  period  wanders 
toward  the  Moon  until  being  thrown  out  of  cislunar  space  (Figure 
30).  The  satellite  travels  50,000  km  in  the  first  lunar  synodic 
month.  When  the  Moon  is  locked  to  the  x-axis  to  eliminate  drift 
the  orbit  is  found  to  be  more  stable,  only  increasing  its  orbit 
size  101  over  the  first  25  months.  But  the  orbit  shape  is  much 
different  making  a  small  3/4  orbit  in  half  of  a  synodic  month 
and  a  larger  3/4  orbit  in  the  next  half  of  a  synodic  orbit. 

This  makes  a  very  irregular  orbit  shape  that  would  not  lend 
itself  to  an  elementary  rendevous  with  a  satellite  in  this  tight 
orbit.  The  smaller  3/4  orbit  has  a  semimajor  axis  of  less  than 
4000  km.  It  is  over  seven  years  before  the  orbit  of  the  satellite 
drifts  out  of  the  quadrant  and  is  thrown  out  of  the  system  shortly 
thereafter.  Note  that  the  orbit  is  unstable  in  VRFB  system.  The 
initial  lunar  conditions  computed  by  B,  M  and  S  are  very  close 
to  those  calculated  by  Kolenkiewicz  and  Carpenter. 
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Fig  30.  Truth  Model  of  BM&S  Orbit, 

One  Orbit  Period 
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With  the  four-body  truth  mode] ,  none  of  the  orbits  found  are 
periodic  bjt  some  do  nearly  close  on  themselves  and  maintain  a 
period  of  one  lunar  synodic  month.  Wheeler's  orbit  is  stable 
enough  to  not  require  a  control  system  to  maintain  it  in  the 
orbit  for  at  least  20  years  (figures  31,  32  and  33).  It  will 
stay  in  the  vicinity  of  L4  for  over  a  century  and  may  never  leave. 
The  transients  are  small  enough  for  space  colonies  but  may  be  large 
enough  to  affect  pointing  in  a  particular  direction  for  military 
applications.  This  would  have  to  be  investigated;  if  the  minor 
perturbations  could  not  be  accurately  predicted  or  there  are 
fast  variations,  then  a  small  controller  might  be  needed  to 
maintain  an  exact  orbit  but  it  certainly  would  not  have  to  be 
large  in  size  or  expensive  in  fuel  consumption. 

Kolenkiewicz  Orbit  I  was  close  to  Wheeler's  but  a  one  percent 
error  in  initial  conditions  is  enough  to  cause  it  to  be  marginally 
stable  (Figure  34).  The  reason  they  didn't  derive  the  same  initial 
conditions  as  Wheeler  is  probably  due  to  the  solar  orbit  used 
which  is  unrealistically  perturbed  by  the  Earth  and  Moon  alone. 

The  same  holds  for  Orbit  II,  which  is  also  marginally  stable 
(Figure  35) .  It  the  initial  conditions  were  derived  using  a  more 
accurate  model,  a  Wheeler  orbit  beginning  on  the  opposite  side 
of  the  orbit  might  be  found, 

Barkham,  Modi,  and  Soudack’s  orbit  is  marginally  stable  and 
not  suitable  for  military  or  civilian  use.  It  is  started  too 
close  to  L4  and  it  has  been  provuid  by  Schechter  and  many  other 
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noted  individuals  that  any  particle  near  L4  will  slowly  be  ejected 
from  the  vicinity  (Figures  3 6f  37  and  38)  ,  Both  of  the  numerical 
orbits  derived  to  check  B,  M  and  S's  orbit  wore  not  truth  model 
checks.  The  first  was  from  the  four-body  perturbation  equations 
but  they  were  simplified  by  driving  the  distance  to  the  Sun  to 
infinity.  The  other  numerical  solution  was  using  a  set  of  three 
body  equations  formed  in  the  same  manner  as  the  four-body  which 
came  up  with  the  loss  than  stable  B,  M  and  S  orbit.  A  stability 
analysis  should  have  been  performed. 

The  fact  that  Wheeler's  orbit  was  stable  in  the  very  restricted 
four-body  model  (circular  Sun  -  Earth  -  Moon  orbits)  and  K  &  C 
and  B,  M  &  S's  orbits  were  unstable  gives  evidence  that  the 
VRFB  model  is  a  valid  simplified  model  to  test  or  search  for 
orbits  which  would  be  relatively  stable  in  the  real  world. 

Wheeler's  orbit  tends  to  support  this.  The  truth  model  also 
decreases  the  VRFB  orbit  by  25%.  Wheeler's  orbit  now  has  a 
semimajor  axis  of  60,000  mi  and  semiminor  axis  of  30,000  mi. 

Figure  39  shows  the  movement  of  one  £>articular  x^oint  on  each 
orbit  for  a  period  of  ten  years  showing  that  the  orbit  moves  very 
little  in  the  quadrant. 

It  was  also  found  during  the  analysis  of  forces  on  a  satellite 
about  L4  that  the  phase  between  the  Moon  and  Sun  is  the  most 
important  factor  in  the  stability  of  the  satellite's  orbit.  The 
Moon  having  its  orbit  perturbed  by  the  Sun  is  the  driving  force. 
Remember  that  Wheeler's  orbit  wi,ll  exist  only  if  the  satellite 
is  initially  injected  into  the  orbit  with  the  Moon  and  Sun  in 
phase  on  the  x-axis. 
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Fig  37 •  No  Drift  Display  of  BM&S  Truth  Model 

Orbit,  3  Orbits 
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Fig  38 •  No  Drift  Display  of  BK4S  Truth  Model 

Orbit,  20  Orbits 
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The  truth  model  itself  could  be  improved  in  accuracy  using 
Ephemcris  Tables  for  accurate  lunar  and  solar  initial  conditions, 
using  the  most  accurate  constants  to  date  from  the  various  authors, 
and  using  more  terms  t.o  higher  powers  for  the  equation  of  center 
which  determines  the  solar  coordinates. 

The  next  step  would  be  to  investigate  out-of-plane  considera¬ 
tions  and  apply  a  third  dimension  to  the  truth  model.  If  the  small 
transients  need  to  be  kept  to  a  minimum  for  the  nearly  stable 
Wheeler  orbit,  then  an  optimal  plane  of  rotation  for  the  orbit 
needs  to  be  found.  It  could  be-in  phase  with  the  plane  of  the 
Moon's  or  with  the  the  plane  of  the  Sun  or  somewhere  in  between. 

The  180°  phase  orbit  could  be  investigated  using  a  Monte  Carlo 
analysis  in  a  very  restricted  four-body  model  to  see  if  a  stable 
orbit  exists  and  then  apply  it  to  the  truth  model. 

Last  of  all,  look  at  attitute  problems  in  the  Wheeler  orbit 
due  to  the  small  perturbations  of  the  nearly  stable  orbit.  The 
initial  conditions  could  be  varied  slightly  and  expanded  to  see 
if  there  is  a  small  family  of  orbits  that  exhibit  stability  in 


that  region. 


Bibliography 


1.  Barkham,  P.G.D.,  et  al^.  "Asymptotic  Solutions  In  the  Many-Body 
Problem"  Part  I  and  II,  Boston:  D,  Reidel  Publishing  Company, 
1975. 

2.  Brouwer,  D.  "Theory  and  Tables  of  the  Motion  of  (588)  Achilles" 
Transactions  of  the  Yale  Observatory,  6  (Pt.  VII) :177-88  (1933) 

3.  Brouwer,  D.  and  Clemence,  C.M.,  Methods  of  Celestial  Mechanics, 
New  York:  Academic  Press,  1961. 

4.  Brown,  E.W.  "On  a  New  Family  of  Periodic  Orbits  in  the  Problem 
of  Three  Bodies",  Monthly  Notices  of  the  Royal  Astronomical 
Society,  71:438-54  (1911). 

5.  Brown,  E.W.  "On  the  Oscillating  Orbits  About  the  Triangular 
Libration  Points  in  the  Problem  of  Three  Bodies",  Monthly 
Notices  of  the  Royal  Astronomical  Society,  71:492-502  (1911), 

6.  Brown,  E.W.  "The  General  Orbits  of  the  Asteroids  of  the  Trojan 
Group",  Astronomical  Journal,  35:69-80,  124  (Aug  1923), 

7.  Brown,  E.W.  "Theory  of  the  Trojan  Group  Asteroids",  Trans¬ 
actions  of  the  Yale  Ob*servatory «  3:1-47,  87-133  (1926). 

8.  Brown,  E.W.  and  Shook  C.A.  Planetary  Theory,  Cambridge: 
University  Press,  1933. 

9.  Charlier,  Meddelanden  Fran  Lunds  Astronomika  Observatoricum, 

No.  18.  . . . 

10.  Columbo,  G.  "The  Stabilization  of  an  Artifical  Satellite  at 
the  Inferior  Conjunction  Point  of  the  Earth  -  Moon  System", 
Smithsonian  Institution,  SAO  special  Report,  No.  80, 

Nov  1,  1961. 

11.  Columbo,  G.  et  al.  "On  the  Librations  of  a  Particle  Near  the 
Triangular  Point  in  the  Semi-Restricted  Three  Body  Problem," 
Smithsonian  Institution,  SAO  Special  Report,  No.  108, 

Nov  20,  1962. 

12.  Cronin,  J.P. ,  et  al.  "Periodic  Motion  Near  Three-Body 
Libration  Points",  Icarus,  3:423  (1964). 

13.  Danby,  J.M.A.  "Influence  of  the  Sun  on  Motion  Close  to  the 
Triangular  Points",  Astronomical  Journal,  70:181  (1965). 

14.  deVries,  J.P.  and  Pauson  W.M.  "A  Study  of  Motion  Near  Libration 
Points  in  the  Earth  -  Moon  System",  General  Electric  Company, 
TIS  R62SD5  (Jan  1962) . 


89 


15.  deVries,  J.P,  "The  Sun's  Perturbing  Effect  on  Motion  Near  a 
Triangular  Lagrange  Point",  XIII  international  Astronautical 
Congress,  Sep  1962, 

16.  deVries,  J.P,  "Motion  of  a  Particle  in  the  Vicinity  of  a 
Triangular  Libration  Point  in  the  Earth  -  Moon  System" , 

General  Electric  Company,  Space  Sciences  Lab  R63SD99, 

17.  Eckert,  W.J.  "The  General  Orbit  of  Hector",  Transactions  of 
the  Yale  observatory,  6  (Pt  VI) ; 161-72  (1933). 

18.  Ellis,  G.A.  and  Diana  A.C.  "Motion  of  a  Satellite  About  the 
Earth  -.Moon  Sextile  Points",  Rome  Air  Development  Commands, 
Technical  Memorandum  RCLS-TM-60-1  (Sep  1960) . 

19.  Feldt,  W.T.  and  Shulman  Y.  AIAA  Journal,  4:1501  (1966). 

20.  Heppendeimer,  T.A.  "Steps  Toward  Space  Colonization:  Colony 
Location  and  Transfer  Trajectories."  Journal  of  Spacecraft 
and  Rockets,  15:305-12  (1975). 

21.  Herz,  H.G.  "On  the  Theory  of  Trojan  Asteroids",  Astronomical 
Journal,  50:121-25  (1943). 

22.  Huang,  S.S.  "Very  Restricted  Four-Body  Problem",  NASA,  TN 
D-501,  1960. 

23.  Kamel,  A. A.  "Perturbation  Theory  Based  on  Lie  Transforms  and 
its  Application  to  the  Stability  of  Motion  Near  Sun-Perturbed 
Earth  -  Moon  Triangular  Libration  Points",  NASA,  CR-1622  (1970). 

24.  Katz,  J.I.  "Numerical  Orbits  Near  the  Triangular  Lunar 
Libration  Points."  Incarus ,  25:356-59  (1975). 

25.  Klemperer  W.B.  and  Benedikt  E.T.  "Selenoid  Satellites", 
Astronautica  Acta,  4:25-30. 

26.  Kolenkiewicz,  R.  and  Carpenter  L,  "Stable  Periodic  Orbits 
About  the  Sun-Perturbed  Earth  -  Moon  Triangular  Points", 

AIAA  Journal,  6:1302  (1968). 

27.  LaGrange,  Oeuvres  Completes,  Paris:  Gauthiers  -  Villars, 

Vol  VI: 277  (1867) . 

28.  Linders,  L.J.  Arkiv  for  Mat,  Astronom.,  OCH  FYS.,  Stockholm: 

SO.  VET.  AK.,  Vol.  4,  No.  20  (1908). 

29.  Michael,  W.H. ,  Jr.  "Considerations  of  the  Motion  of  a  Small 
Body  in  the  Vicinity  of  the  Stable  Libration  Points  of  the 
Earth  -  Moon  System",  NASA,  TR  R-160,  1963, 


90 


w 


30.  Musen,  P,  and  Carpenter  L.  "On  the  General  Planetary 
Perturbations  in  Rectangular  Coordinates",  Journal  of 
Geophysical  Research,  68:2727-34  (1963). 

31.  Plummer,  H.C.  "On  Periodic  Orbits  in  the  Neighborhood  of 
Centers  of  Libration,"  Monthly  Notices  of  the  Royal 
Astronomical  Society,  62:6-17  (Nov  1901). 

32.  Rabe,  E,  "Determination  and  Survey  of  Periodic  Trojan  Orbits 
in  the  Restricted  Problem  of  Three  Bodies",  Astronomical 
Journal,  66  (No.  9) : 500-13  (Nov  1961). 

33.  Rabe,  E.  "Additional  Periodic  Trojan  Orbits  and  Further 
Studies  of  their  Stability  Features",  Astronomical  Journal, 
67:389-90  (Aug  1962), 

34.  Rabe,  E.  and  Schanzle  A.  "Periodic  Librations  About  the 
Triangular  Solutions  of  the  Restricted  Earth  -  Moon  Problem 
and  their  Orbital  Stabilities",  Astronomical  Journal, 

67  (No.  10) :732-9  (Dec  1962). 

35.  Roosen,  R.G.  "Objects  in  the  Vicinity  of  L4  and  L5", 

Icarus ,  9:429  (1968) . 

i 

36.  Schechter,  H.B.  "Three-Dimensional  Nonlinear  Stability  Analysis 
of  the  Sun  -  Perturbed  Earth  -  Moon  Equilateral  Points", 

AIAA,  Paper  67-566  (1967) . 

37.  Schutz,  B.E.  "Numerical  Studies  of  Solar  Influenced  Particle 
Near  L4",  Unpublished  Dissertation.  University  of  Texas. 

38.  Siferd,  R.E.  "Generating  Periodic  Orbits",  Unpublished  Thesis. 
Wright-Patterson  Air  Force  Base,  Ohio:  Air  Force  Institute 
Of  Technology,  December  1965. 

39.  Smart,  W.  "Libration  of  the  Trojan  Planets",  Paper  I,  Memoirs 
of  the  Royal  Astronomical  Society,  62:79-112  (1918). 

40.  Steg,  L.  and  deVries,  J.P.  "Survey  of  Motion  Near  Libration 
Points",  Space  Science  Review,  5:210  (1966)  . 

41.  Stumpff ,  K.  "On  LaGrange's  Theory  of  the  Three-Body  Problem", 
NASA,  TN  D-1417  (Jan  1963) . 

42.  Szebehely,  V.G.  Theory  of  Orbits,  New  York:  Academic  Press, 
1967, 

43.  Tapley,  B.D.  and  Lewallen,  J.M.  AIAA  Journal,  Vol.  2  No.  4 

(1963) .  '  “ 

44.  Tapley,  B.D.  and  Schultz,  B.E.  AIAA  Journal,  Vol.  3  No.  10 
(1965)  . 


I 

3 


u 


91 


45,  Tapley,  B,D,  and  Schutz,  B.E,  AIAA  Journal,  6:1405  (1968). 


46,  Thuring,  B.  "The  Trojan  Libration  Periods  of  Orbits  of  Various 
Amplitudes",  Astronomische  Nachrichten,  243:181  (1931), 

47,  Thuring,  B,  "Seeking  Solutions  to  the  Theory  of  Movement  of 

the  Sun  -  Jupiter  Group",  Astronomische  Nachrichten ,  279:217- 
30  (Sep  1950)  ,  ' 

48,  Thuring,  B,  "Program  Formation  of  Libration  Station  Keeping", 
Astronomische  Nachrichten,  285:71-80  (1959). 

49,  Wheeler*  J.E,  "Determination  of  Periodic  Orbits  and  their 
Stability  in  the  Very  Restricted  Four-Body  Model  in  the 
Vicinity  of  the  Lagrangian  Points  L4  and  L5",  Unpublished 
Thesis.  Dayton:  Air  Force  Institute  of  Technology,  Dec  1978. 

50,  Wilkens,  A.  "Methoden  Zur  Ernittlung  Der  Speziellen  Und  Dei 
Absoluten  Koordinatenstorungen  Der  Planeten  Der  Jupitergruppe" , 
Astronomische  Nachrichten,  206:149-56  (Apr  1918). 

51,  Wilkens,  A.  "Integration  of  the  Differential  Equations  of  the 

Jupiter  Problem,"  Astronomische  Nachrichten,  206:149-56  (Apr, 
1918) .  , 

52,  Wilkens,  A.  "Concerning  the  Libration  Theory  in  the  Sun  System", 
Astronomische  Nachrichten,  227:337-48  (May  1926). 

53,  Wilkens,  A.  "Understanding  of  the  Sun  -  Jupiter  Theory", 
Astronomische  Nachrichten,  245:229-64  (May  1932). 

54,  Willard,  H.R.  "On  a  Family  of  Oscillating  Orbits  of  Short 
Periods",  Monthly  Notices  of  the  Royal  Astronomical  Society, 
73-471-84  (1913) . 

55,  Wolaver,  L.E.  "Effects  of  Initial  Configurations  on  Libration 
Point  Motion."  Progress  in  Astronautics,  17:75:99  (1966). 

56,  Wolff,  C.  et  al.  "Cloudlike  Objects  at  Triangular  Points", 
Science,  157:427  (1967). 


k 


92 


APPENDIX  A 


93 


Comparing  equation  (2)  to  equation  (1)  with  the  assumptions 


that  the  mass  of  the  ith  body  remains  constant  and  drag  and  other 


external  forces  are  not  present,  we  obtain 

El  =  -  G  ;  <£,i) 

3=1  r  3 

3*1  ji 


(3) 


If  m2  is  a  satellite  and  m^  is  a  planet,  writing  equation  (3) 


for  both 


r.  =  -  G  Z  _j  (r . . ) 

_1  J-2  3 

r . . 


n 


m. 


r2  =  -  C  Z  _j.  (r  ) 

3=1  _  3 

j*2  rj2 


(4) 


(5) 


And  since  r12  =  r  2  -  r  ^  so  that  f  =  r  2  "  Substituting 

* 

equations  (4)  and  (5)  into  this.  last  equation 


fo.2 


=  -  G  E  ^  (^  )  +  G  E  _j  (r  ) 

3-1  3  ^  3=2  r  3  jl 

j^2  rj2  rjl 


Expanding 


-12  ~  ' 


(— 12^  +  G  ^  V 

r  3  3=3  3 

r12  rj2 


n 

-  G  m_  ,  .  G  E  m.  ,  . 

—1  (^1}  "  j=3  -I  (^1> 

r  3  r  3 

r21  rjl 


Combining  the  first  terms  in  each  bracket  since  r^2  =  -  i: 


ri2  =  -  G  (mL  +  m2)  „  s  G  ■  -j  /  ^2  . 


3=3 


’12 


r  3  r  3 

Cj2  rjl 
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ttll 


=  f-  G  (m.,  +  m0)  .  .  _  G  jn.  /r  -  r .  r.  ^-r. 

12 _ i— -i-  (£12)  -  S  3  /=i - 3L  _  ^__ZL 


'12 


3-3 


"2  j 


Consider  an  Earth  -  Moon  -  Sun  system  with  a  colony  satellite 
in  the  vicinity  of  the  lagrangian  point  L4  or  L5.  Refer  to  Figure 
1.  Assume  the  following: 

1)  Coplanar  Orbits 

2)  Mass  of  the  satellite  is  negligible  compared  to  the 
masses  of  the  other  bodies. 

3)  The  Sun,  Earth,  and  Moon  are  considered  point  masses. 

4)  Ignore  effects  of  other  planets,  upon  the  satellite. 

From  Figure  1 


r  = 

x  a. 

+  y 

a_ 

w 

— s 

s  -1 

s 

—2 

r  = 

o 

—e 

r  = 

x  a. 

+  y 

a„ 

— m 

m  —1 

m 

—2 

r  = 

x  a. 

+  y 

a„ 

— c 

c  —1 

c 

—2 

r 

=  (x  - 

x  ) 

a. 

+ 

(y 

-  y  ) 

a_ 

— m s 

m 

s 

—1 

m 

s 

—2 

r 

=  (x  - 

X  ) 

a. 

+ 

(y 

-  y  ) 

a. 

— cs 

c 

s 

—1 

c 

s 

—2 

r 

=  (x  - 

X  ) 

aT 

+ 

(y 

-  y  ) 

a0 

—cm 

c 

m 

—1 

c 

m 

—2 

and 

r  2 

=  X  2  + 

y_2 

r2  =  x2  +  y2 
s  s  s 


r2  =  x2  +  y2 
c  c  c 


r  4  =  (x  -  x  )  *  +  (y  -  y  ) 
ms  ms  Jrm  a 
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l 

! 

I 


r  •i=(x  -  x  )  ^  +  (y  -  y) 
cs  c  s  c  s 


r  ^  =  (x  -  x)z+(y  -  y  )  ‘ 
cm  cm  cm 


Using  the  perturbative  function  equation  (6)  for  lunar  motion 


1  «*  H 

where  r.  0  =  r 
— 1/  in 


„  G(m+m)  _  / r  -  r 

r  =  -  e  m  r  -  G  m  /  — m  — s 


Where  constants  m  =  p 
m 


m  =  1  -  p 
e 


G  =  1 


So  equations  for  lunar  motion 


««  X  /XXX 

X-D1R:  x  +  m=  -m  /  m-s+s 
m  —  si  - -  — 

*r  3  l  r  3  r 

m  \  ms  s 


Y-DIR:  y  +  y  =  -  m  /y  -  y  +  x  ' 
m  m  s  /  _m _ s 

r  3  1  r  3  r  3 

m  \  ms  s 


Using  equation  (6)  for  the  colony  satellite  motion 


r=-G(m+m)  Gm  /r  -  r 

-c  e  c  r  -  S/-C-E 


-x"- 


-Gm  /  r  -  r  ,  r 

m  /  ~m  -  — t 


Using  the  same  constants  and  m^  negligible  compared  to 


L  +  U  -  n>  ^  ,  _  »s  /r. 

r  3  l  r  3  r  3) 

c  V  cs  s  / 


r  -  r 
— c  — m 


X-DIR:  x  +  l-  p  x  =  -  m  /x  -  x  x  \  /x  -  x 
c - c  s/c  s+s\-Vvc  m 

~  r>) 

\  cs  s  /  \  cm 
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(9) 


y-dir;  y_  +  l  -  u  ..  _  _  /y_  -  y.  .  y_\  /  y_  -  y_  .  y_N 

c  — — — y  —  m  /  c  s  +  s  \  .  /  c  m  +  m 

ro3  c  r>r  ~  73 

\  cs  s  /  \  era  m  . 


Cast  these  four  second  order  differential  equations  (7f  8 
9  and  10)  into  eight  first  order  state  equations 


So  that  the  following  differential  equations  result 


x,  x 

1  =  m 


2  * 


x,  x 

3  =  m 


4  = 


5 


6 


7  = 


'8  = 


In  tejrms  pf  the  state  variables,  the  equations  of  motion 


become 


k  =  -  (1  -  u) 

3 

r 

c 


(12) 

(13) 


(18) 


(19) 


Where 

3  _  .  2  2.  3/2 

r  =  (x  +  x  ) 
c  5  6 


r  3  *  r(x  -  x  )2  +  (x  -  y  ) 2]  3/2 

CS  5S  OS 

3  r/  .2  .2,  3/2 

rcm  =  [(X5  '  V  +  {X6  "  X2)  ] 

r  3  *  (x  2  +  y  2,  3/2 

a  s  s 

3  .  2  ^  2,  3/2 

r  =  (x  T  +  x  ) 
in  x  z 


3  r,  >2  ,  ,2,  3/2 

rms  *  t(xl  "  V  +  x2  "  V  3 


(20) 

(21) 

(22) 

(23) 

(24) 

(25) 
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APPENDIX  B 


DERIVATION  OF  THE  EQUATION  OF  CENTER 
Sun  coordinates  are  needed  for  the  satellite  equations  of 
motion  in  the  truth  model.  Motion  of  the  Sun  is  defined  by  the 
equation  of  center. 

A  fixed  system  of  rectangular  coordinates  with  the  origin  in 
the  center  of  the  Sun  has  equations  (Ref  3)  of  motion 


d2x  ,  x  3R 

7T  +  11  “3  =  w 

dt  r 


(1) 


+  u  M 


dt 


3y 


d2z  ,  z  3R 


where  m  =  1,  ra  is  the  disturbed  planet,  m'  is  the  disturbing 
s 


planet,  y  =  G  (1  +  m)  and 

•  2  '  2  ' 

R  =  [(x  -  x)  +  (y  -  y)  +  (z  -  z)  ] 


I  f  I 

21  3/2  x  x+y  y+z  z 
,3 
r  ' 


Perturbed  coordinates  are 


x  =  x  +5 

O  X 


x  =  x  +6. 
O  X 


y  =  y  +6 
o  y 


y  =  y  +5. 
o  y 


z  =  5 


z  =  6. 


with  the  x-axis  directed  toward  the  perihelion.  Substracting  the 

equations  of  elliptic  motion  from  (1)  yields 

2  x 

— ,  (<5x)  +  -  (x  +  5x)  -  y 

dt2  r3  0  r  3  6x 

o 

(same  for  y  and  x  components) 
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Develop  the  second  item  in  the  left-hand  members  in  powers  of 

$  .  6  ,6 

x  y  z 

2 

Ex:  — 2  (Sx)  +  ^  3  (6x)  -  ■  - -y-  (xq  5  x  +  6  y)  =  Gx 

dt  r  r  A 

o  o 

and  expand  the  partical  derivatives  of  the  disturbing  function  in 
Taylor’s  series  in  terms  of  R 


2  2  2 
Ex:  G  =  3r  ,  3  Ro  (  ,  3  Ro  ,  .  3  Ro 

x  +  -  ox  +  - - ; -  o  y  +  - - - —  6z 

3x  »  3x  dy  3x  3z 

o  ax  o  o  o  o 

o 


+  3  Ro  ,  ,  ,  3  Ro  „  ,  ,  a  Ro  ,  . 

3x  3x  3x  oy  ’  1  3x  3z  ' 

o  o  o  o  o  o 


+  H  [(9  Xo  - 

1^5  o  )  <5  x  +  (3  yo 

2  5 

y 

2  7  : 

L 

O 

ro  ro 

+  (3  Xo  -  15 

Xo  yo  )  6  y2  +  3  Xo 

2  5  2 

7  2  5 

r 

o 

r  it 

o  o 

Integrating  the  expanded  equations  (first  order  approximation  to 

solve  second  order  of  disturbing  forces,  etc.)  to  find  a  formal, 

explicit  solution  yields  expressions  containing  the  terms 

6x  =  C  +  C  t  +  C,  cos  L  +  C_  cos  2  L  +  C,  cos  3  L  +  . . . 
x  ox  lx  2x  3x 

+  S^x  sin  L  +  S2x  Sin  2  L  +  S^x  Sin  3L  +  ... 


(6y  is  similar) 

2  3 

where  C  and  s  are  numbers,  secular  and  terms  factored  by  t  or  t 
are  ignored.  Time  is  the  independent  variable.  Converting  to  a 


true  anomoly 


6v  =  --  (x  5y  -  y  5x) 

2  O  Q 

r 

o 


=  C  +  C  t  +  C,  cos  L  +  C„  cos  2  L  +  C,  cos  31+,.. 
o  I  2  3 

+  sin  L  +  S2  sin  2  L  +  sin  3  L  +  ... 

Expressing  a  finite  increment  of  v  as  a  function  or  arbitrary 
increments  of  the  four  elements  u>,  n,  e,  S3  and  -  u  =  L  =  nt 

for  the  Sun  as  a  disturbing  body 

v  =  nt  +  (2  e  +  7  e3)  sin  nt  +  (7  e2  -  ~  e4)  sin  2  nt 
4  4  24 

+  y|-  e3  sin  3  nt  for  3rd  order  eccentricity 

Use  this  to  calculate  the  Sun  coordinates  starting  at  perihelion. 

6x  *  r  cos  v  <5y  =  r  sin  v 
o  o 


therefore 


where 


X 

=  y 

X 

+ 

r 

cos 

V 

s 

m 

s 

=  y 

y 

+ 

r 

sin 

V 

S 

m 

s 

2 

r  =  A  (1-e  )/(l  +  e  cos  v) 

s  s  s  s 
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APPENDIX  C 


COMPUTER  PROGRAM  TO  CALCULATE  ORDITS 
IN  THE  RESTRICTED  FOUR-BODY  TRUTH  MODEL 
This  appendix  contains  the  computer  routines  used  applying 
Wheeler's  and  Kolenkiewicz  and  Carpenter's  orbital  data  to  the 
truth  model.  Barkham,  Modi,  and  Soudack's  routine  is  not  displayed 
because  of  its  similarity  to  Wheeler's  routine.  The  computer 
language  utilized  was  fortran  Extended  Version  IV,  and  all  work 
was  accomplished  on  the  AFIT  CDC  6613  and  CYBER  74  computers , 
Several  comment  statements  have  been  employed  to  aid  the  reader 
and  smooth  the  flow  of  the  program.  ODE  was  the  integration 
package  utilized  in  the  program.  For  the  integration  steps,  one 
time  unit  corresponds  to  4* 3483774  days.  Data  cards  were  input 
into  the  Kolenkiewicz  and  Carpenter  routine  containing  the  a  and 
6  coefficients  from  Tables  3  and  4. 
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PROGRAM  GRAOU  74/7*.  OPT=i 

Wheeler  Truth  Model  with  Subroutines 


PTN  **•  7 


5r  rruxn  ivioaei  wiun  oumuuoinca 

PROG'-  A  I-  GRAOU  (INPUT  =  /A  '  ,  OUTPUT  ,TA  F  E>  =  INPUT  ,  T  AFF6  =  OUTPUT  , 
EXTERNAL  Fl 

ki  O  AC  .  >  *  .  411.  CT 


°L  -  “ 


COMMON  MS  »AS»  ;S»f:  j»  1U»FI 
REAL  MS»NS»MU 

DIMENSION  X(8),XM(1300) ,YH(l3Ci>, VXN(13Cf> ,VYM( 
-YC  (13  to  ,  VXC<130;>  ,  >/YC(13&:>  f  TMUJU) 

DIMENSION  IWO*K(r  j-1)  .WORK^dO) 

DIMENSION  XW(  13  If  )  »YW(133u)  »  VXW  ( 1 33  J)  ,VYW(  13l0) 
0IMENS10N  AA (4? 1)  *93  4:  0) 

MS =32 6901* 12 


y/YM(l?Uu)  ,X0(1T 


\jx  n  o  jl  wiy  ««  it-  it  j 

MS =32 690. 1.12 
AS  =30  8 . 620284  9 
ES  =  .l  166 
NS  =  .l  7-46-13 
MU= • t 12139609  ♦ 

PI  =3.  141*926936 
XMIN= It « 

YM1N= 1C. 

XMAX=C. 

YM  AX= 9. 

X ( 1)  =  -  .9S22C573479 

X(2)=C 

X( 3 ) = C 


X( 3 ) = C 

X(*»)  =  -l»ulG56l6  13':3fi 
X(E)  =  -. 7363274299° 

X(6>= .61io9639Gd9  * 

X(7)  =  -X(f.  )♦.  V94(  i3194°/i. 1606489351 
X(  )  =  X  (3)  4.224  33  i  0."  7  58/1.0836489351 
N£QN=  6 


T=C 

TOUT=2.u-  PI*  1.0  5:  548  93F 1/13  3. 3 
OELT= 1  OUT 

RELERRs.^CCCO  JOOCl 

A3SEF  Rs.SOUOJ  l^PCl 

IF  LAG  =1 

KCUNT  =1 

Trt  (KCUNT)  =T 

Xrl  (KC  UNT)  =X<1) 

YM (KCUNT) *X(2) 

VXMtKOUNT )=X(3) 

VY  M (KOUNT ) =X (4) 

XC (KOUNT) =X(9) 

YC (KCUNT) =X(G) 

VXC (KOUNT ) =X  (  7 ) 

VY C  (KOUNT  )  =X  ( 6 ) 

A*ATAN(X(2)/X(1)) 

IF  (X(  D.GT.L.1)  A  =  A  +  FI 

IF  (X(  1)  .LT.O  .  3  .  AND.  X  (2)  .  GT.o  .  0)  A  =  A*2.-*PI 

FS=NS*T*(2*ES-. 2  *ES ** 3) *  SIN (NS* T ) +1. 26*  ES  *  2*SIN(2*NS*T) 4(1 
-♦tS** 3*SIN<3*NS*T) 

FS*FS  +PI 

1  .  if  4  ft  a.  cci  r  nc  i  cz  \  \ 


FS*FS  +PI 

RSP*AS* ( 1. u-E5*  *2 ) / ( 1. p  *ES*C0S (F» ) ) 
XS=hU4  X(l) +RSP*  COS (FS) 

YS=MIJ*X(2 )  4RSP*  SIM(PS) 

PRINT #  t**T  =  -,f,"  XS  -  "tXS,'*  YS*  ",YS 
XH  (KOUNT)  =X(>)*  COS  (A  )4V  (6)*SIN  ( A)  +M'J 
YM (KOUNT) *-X(5) •$ IS ( A) 4X  (6>*C0S(*> 
VXW(KCUNT)  =  (X(7) +YW( KOUNT) ) *  1 • 0 33  8h3 3 31 


51 


53 


/d 


7  5 


33 


35 
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VY  W  (KQUNT  )  =  (X  ( 3  )  -XW  (  KOUKT  )  )*  1. 4  83  3^35  331 
3  CALL  GCE(F1, NEON, X, T  ,TC  UT,RELERR,  ATSERR.IFlAG, WORK, IWORK) 

lF(IFLAG.EQ.l)  GO  TO  F 

IF  (IF  LAG. ME. 2)  GO  TO  2 

KOUN’  =KOUNT  +  l 
IF  (KCUNT.E9.25n  GO  TO  4 
TO U7  =  TOUT ♦DELI 
GO  TO  1 
t  continue 

PRINT  *  » "MOON  0%1IT" 

00  6  1  =  1*  KOUNT 

WRIT? (6,103)  TM(I),XH(I),YM(I),V<M(I) ,VYM(I> 

19C  FORMAT  (13X,"T=  *•,  F14  .9  ,HX,"X=  " ,  l  IS.  1 1 , 5  X,  “Y*  ",  G16.ll  ,5X,"VX 
-Gl5.11,5X,“VY  =  ", G13.il) 
t>  CONTINUE 

PRINT*  ."SATELLITE  OR  FIT" 

DO  7  1=1. KOUNT 

WRiTE  ( l ,  3G 1)  TM  (I )  ,XC(T)  ,YC(I)  ,V<C(I)  , VYC ( I) 

30b  FORMAT (1«X,"T=  ", F14 .9 ,4X  , "X  =  " , 1 13. 1 1 , 3 X, "Y=  ",  G13.ll ,5X, "VX 
-Glo.ll,5X,"VY=  ", G18.il) 

7  CONTINUE 

PRINT*, “WHEELER  OFT  IT" 

00  9  l=i,KJUNr 

WRITE  (6,3 9 2)  TM(I),XW(T),YW(I),VXW(X) ,VYW(I) 

500  FO  RMA  T ( 1 L  X ,  "T  =  ", F14 .9 ,4X  ,"X  =  ",G18.ii,5X,"Y=  ",G16.11,5X,"VX 

-Glc.l 1,5X,"VY=  ", G13.il) 

9  CONTINUE 

N=lil 

Do  12  1=1,  N 
AA (I) =XW ( I) 

BS  (I)  =  YW  ( I) 

12  CONTINUE 

MM=N 

CALL  FLOT  SC ( A  A , 9F ,MM ) 

2  PRINT  * ,"I FLAG=  ", IFL  AG 

3  CONTINUE 
STOP 
ENO 


SUBROUTINE  PLOTXY  74/74  UPT=1 


FT  N 


1  SUBROUTINE  PLOTXY(X,  Y, N , XMIN , XMA< , Y U N, YMA X, 1ST ) 

DIMENSION  X  (N)  ,  Y(N) 

CALL  PLOT  O.Cy-l.  0,-3) 

CALL  PLOT (j.C,0.i,-3) 

5  XFAC=1L.0/(XMAX-XMIN  > 

VFAC=1C.  j/(YMAX-YMIN> 

00  959  1=1, N 
XX*(X  ( I)  -XMIN)  *  XF AO 
YY  =  (Y  (I)-YMIN)*  YFAC 

1J  CALL  SYMBOL(XX,  YY,T. 07f , 3,0.u,-1* 

999  CONTINUE 

CALL  PLOT (10., 0 .0, -3) 

IF  (IST.EQ.O)  RETURN 
JJ*0 

IS  CALL  PLOT E(JJ) 

RETURN 
ENO 
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SU?*OUTI.NE  Fi 


7*/7h  0°T  =  i 


FT  N  4*  7*<m  t> 


1 

5 

13 

15 

23 

25 

33 

SUR ROUTINE 

1 

5 

11 

15 


SU°KCUTInE  Fl(T,X,OX> 

C0HNON  MSf^S»ES,NS,H'J,Pl 
PEAL  PS,NS,rtU 
DIMENSION  X<6)  ,  DX  (3) 

Fb  =NS  ♦!  ♦  <  23*  IS  **I)  *SIN  (NS*T  )  +1. 2.5*  E$  ’  *  2*Z  IN  <  2*  MS4  T>  f  (?3i 

-*£S*<  3*SlN(3*MS*r> 

FS  =  FS*PI 

PSP*AS*  (1  «u-E3»  *2)  /  (  i»7*e S* COS  (PS  ) ) 

XS  =  MU*X(1)  ♦Ro°*lCuS(*S) 

Y$=hU*X(2)+RS»*S..N<PS) 

RMS-i  (XfD-XS)*  *?  +  (*.  <2)-YS)«*2)*>  ,5 
RMMX  <1)**2+X<2)  -*2>  ‘*,5 
FS=(XS**2+YS**2  >*  \3 
RC  =  (X  (t)  *  *  2+X  (£• ) '  '  2)  **.5> 

RCS=(  ( X  ( ?  )  -XS)  *  *2  +  ( X  (E)-YS)4*2>4>.5 
FCh=<  (XCb)-Xti)  >* *  2+  (X(6)-X(2))*»2)  **  ,5 
OX (1) =  X(  3 ) 

OX  <2)  =X(t ) 

OX  <3)  =-X<  i>/Rt-*?-4SJ,<<X(l)-XS>/?HS4<' 3+XS/-Sv*3) 

OX  CM  s-X(2)  /Ri»  *3-M3*{(X(2)-YS)/?M3**3  +  Y£/r<:S#*3) 

CX  (?)  *XCT  ) 

OX  <b>  *X<3) 

DX  (7 )  “X  (3)  /F.r**3-MS»  (  ((X(3)  -rfu*  X  (1)  -RSF4COS(FS)  )  /HCS 

■M  (MU+X(i>fRS°‘COS(FS))/RS-*3))-4U*n  <X(>)  -X(1))/RCM“  3) 
•♦X(1)/FM'*3> 

OX(6)«-(l,j-rtJ)^X(5)  /Rr*43-MS*  <  (<X<3)  -MUM  < 2) - <3P* SIN (FS)  )/RCS 
■♦(  (Ml>*  X<2  )  +RS°*  Sirups)  )  /RS»*3))  — 1 J*  ( (  (X(  fc>)  -X  (2)  )  /RCM**3> 

'♦X  (2)  /RM**3) 

RETURN 

END 


PLOTSC  74/7**  OPT  =1 


FTN 


SUBROUTINE  PL 3T SC (X , Y,M) 

OINENSION  X(M>  ,  Y(  M> 

CALL  PLOTS <3l!) 

CALL  PLOT  (Q.  »-3  *3) 

CALL  PLOT (0., .5,-3) 

CALL  PLOT(1.2?, .;,3) 

CALL  PLOT (1.23,9.5,2) 

CALL  PLOT  (7,23, 9.  ",  2) 

CALL  PLOT (7,25, ,?,2) 

CALL  PLOT (1.25 , .5 ,2) 

CALL  PLOT  <1.-/3, 1.  0,-3) 

CALL  SCALE<X,5.,M,i) 

CALL  $CALE(Y »3.»K,i> 

CALL  AXIS (0., J. ,?HX-AX7S,-6,5.,3» ,X( H  +  l) , X <M+2)  ) 
CALL  AXIS<3.,0.,tHY-AXIS,6,8.,9l.,Y(4Ll),Y(M+2)) 
CALL  L1N£(X,Y,M,1,-1  ,3) 

CALL  PLOT  E(Q) 

RETURN 
END 
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Kolenkiewicz  and  Carpenter  iruth  :<.odel 

Fr\CG  A  f  G^AC'Jd  K'jrs /8' ,  0UTPU1  ,T»c:;sI  »*I 
FX^E'  NAL  Fi 

COMMON  -IS  ,  AS,£S,>  8,MU,M,  AX  (lb)  ,lr'(i;  )  ,  3C{  1  =  )  , 
R'ttL  M«,i,5»MU»u3,  Cl,  mmm 

Gi «=r  siu.,  x ( - ) ,  x*  <i  ^  ,  Yhd  3.  -  ),  i  <  1  i  ) ,  7v-  • 

-yc  (i  ?o( ) ,  vxc<i?  r  > ,  j  YCfi?j«.)  ,ii(i?;j) 

Jir.tr  mon  iwoik  (r  .  :>  ,wr  r > 

CIM?  $10:«  <«(13  1'  ),YXC3j3),VXVI(13il),VY«(:1>  ' 

DIM)  SICi.  A  A  (2i  ~)  ,  ■P  (2<  3) 

11  R{hO(-,*)  *  A"*  »°  0, 

IF  ( £0  F  (v )  .,T“.  J)  GC  TC  3 
CALL  EhXL  £T ( 1 J » 1) 

+*5  =  32  b*:l  2 .42 
AS  =3.'  £-.32  32 
FS  =  .  It 3 
NS=..  743.13 

hU=.v  121‘-  '>.7 

PI-3.  1-1.9261  3c 
AC  =  3.  313  >123.  ’ 

A*  =  3.  i?Hi7'io.  j 

M  =  (i  '<*2.3  V  PI)  /(~Sf  2r.  +  2,.*  2 '  n  .  1<  J  .  •  3-1  . 

M  =  (?9.ii3uG332/  (2*’I)  )•  24.*bl  .  +  6)  . 

N3  =  • ..  i  j  ♦^F'T3cS 

alfh;  =  „<<  1) 

BE  1 A=  PC ( 1 ) 

AO 07= t 

90  07=  t  » 

DO  L  1=1.13 

ALFH/  =ALaHA+A<(I+i)»  C0P  CI*1.  Oh71?78>  UD(3n)  r 
9ETA=  E  ET*'  +3C  (I  +  l)  •  COS(T*  l.u'/in  )  +3?  (i  +  l)  M' 
A007=  ACOT-AK  (I+i>«  I*  (N  3 -Nl)  +  SIN  (I  *  1. 1  +(  1 ~  .* '  >  +  A  ' 
-COSC  *1.  ) 

BD0T=  EC0T-8C  (I  +  l)  *I‘  (N7-N1)  *  SIM(  T  •  1.  [  *71  9  ( ' )  +  ar  (  . 

-cosa  +  i.  ^/ip’?  ) 

9  CONTINUE 

AlFHA=ALPHA*.  J-3  V  31 
5E7  A=  E ETA*  •  u0  Jn  71 
AD0T=  AC0T*.C  j  J£  II*  MS 
BC07=6C0T*.;«j  JC  11J  M 

PlINT'iALFHA,’*  "  » 3E  7A  » "  ",A00T,"  ",00(1 

X(1)=1.G-. 10779420331 
X  ( 2 )  =  C  • 

X  (  3 )=  C  • 

MMU=(  N3+*  2)  *  (AM  “  3) 

X(h)=  (4,;857£3b"-c /’M?+X(l)  >  +  (MMU/AH)  «  +.9*  (NN/(  '-; 
CMU=f  N3+*  2)*  (AC+*3) 

X(3>  =  (1  +  4LPHA)+  (AC/AM)'  C0S(  1.1471  573)  +  3  r  T  A  •  (l./<* 
•  (  ( CMC  /AC)  +*•£)*  (MN/AH)-  ( -SI  N  ( 1.  0. 71  3"  fc )  ) 

X(b)=  (1+ ALPHA)  *  (AC/ A  V)  *SIN(  1 . 0  471  5  7S)  +  3ETA  •  (1./  C 
-((CMU/AC)".E)  *  (  M N7  A  w)  '  COS(  1  •  U47L  973 ) 

X<7>=  (A00T-(0ET  A/N'7)  *  (0 MU/AC*  *3 )  *  NN)  f  (AC/A.-1)* 
-COS(1.C4?1970)  +  (PCOT/(N3*NN)  +1.+ALP3A)  MCMU/AC) 
-(-SIN  (1. 0*71976)) 

X(  7)  =  X  (7) /1.C30 89331 

X(t>  =  ( A0GT-(8ET  A/'H)  *  (f  MU/AC*  +  3 )  *  NS)  *  (AC /AY,)' 

-SIN  (5  •  w47 1976)  +  ( F 00 T  / ( N  3*  NN )  +  i.+UDHA  )  *  (  CHU/AC) 
-COS (1  •  t*/ 197b)  ^08 


X(  c  )  =  X  (3)  /l.f  '3l  1-  PV?'ii 

PkINI  *  ,"X  ♦  ”,X<),"  X5  ",X(5),m  Xu  ",X(f> 

-.1  Pk1kt*,“x7  **,x  (:),•*  ya  “,xm 

MCN=£ 

t  =  >. 

Ti)UT=2.l»'PI*l.i:-1  ^•«r?*l/3C.C 
P£  LT=  TOUT 

vi  ’  tLtl  !t  J  01  .1  j  1  -  1 

A3SE-  R*.  jfln-H  JO  V.  1 
3FLAr  =1 
<0  UNT  =  1 

1  7M  IKljtt.T)  =T 

70  *,t  (KC  UKT)  =  XCi) 

YM  (KC  UNT)  =X(2) 

VX  r (K  L-UN1  ) =X ( 3 ) 

VY  .V  ( Y-C UNT  )  =X(4) 

XC  (KCUNT)  =X  (:•) 

•  -3  YC  (KO  UNT )  =  X  ( *-  > 

VXC ( <  CUNT ) =X ( 7) 

VYC  CUN1  )  =X  (  3 ) 

A=ATA  N  < X (  2)/X(l  )) 

IF  (X(l).LT.b.J)  A  =A  +  DI 

30  IF  <X<  1)  .37  .J.  J.  V'H.X  (2)  .LT.L  .  0)  5=44-2.  -PI 

yu  (<r.uM)  =x(i/)*  crs  ca  )+v (e)»j>iN(A»  +^'j 
YW  (KC  UNT)  =  -X  (  5)  ‘:ri(  A)  4-X  (S)  *COS  (S  ) 

VX w ( K CUNT  )  =  <X<7)  -  X(E>  n  COS(A)  +  (Yf  3)  4-x  (“>  -I'D  ‘STrHA) 

VYW  (KOUN7  )=-(X(X)  Tx  (6)  )“SIN(a)  4-  (X  (3)  +  X  (5)-i3U)*uOS(A) 
c5  5  CALL  CCeitFl, NEON, X,T,tcuT,kELERP,A93EKR,1FLAG, WORK, IWORK) 

IF (IF Lf G. EO.l)  GO  rO  * 

KQ  U  NT  =  K  OU  NT  4- 1 

TF  ( KC  UNT  • tQ. 25u )  GO  TO  '> 

Tu  UT=  1  CU  T  ♦QELT 
30  GO  TO  1 

*>  CO  NT i  N  L  c. 

P:<*Nf  *,"MCON  ORRIT” 

00  fc  1=1,  KOUUT 

WRITE  (f,icj)  TM(J),XMD  ,YKCI) ,VXH (I)  ,VYH<I) 

03  M1  FOPMAT(i.X,"T=  ",  4  .9  ,4  X  ,"X=  ",  3 13 . 1 1 , 5  X,  "Y=  ” ,  G 1 8 •  1 1 ,  i>  X ,  "  VX 

-G16.1 1 ,5X,"VY=  ", G13.ll) 

6  CONTINUE 

PRINT*, “SATELLITE  OR=TT” 
nO  V  I=l,KQUNr 

1.3  WRITE  <t, 300  TW  (I),XCC)  ,YC(I),V<:(I)  ,VYC(I) 

330  FOLMA  T  (1 jX,”T=  ", F14  .0 , 4 X  ,”X  =  ", 313.11, EX, "Y  =  ", GiO. 11 , 5X , "VX 
-Glb.ll,-»X,"VY=  ”, G18.ll) 

7  CONTINUE 

PRINT*, "WHEELER  OnRIT” 

1.3  CO  9  1=1 .  KOUNT 

WRIT?  (t,1;  0  0)  TH  (I  ),XW(T)  ,  YW(I)  ,VXW(I)  ,VYW<1) 

E-  JO  FORMA  T  (1 '  X,”T=  ”,  F14  .9, 4X,"X  =  ”,313.11,  EX,  ”Y=  ”, G18.ll  ,5X,”VX 
-G13.il, 5X, "VY*  “, G19.il) 

9  CONTINUE 

111  N=  111 

00  12  1=1, N 
AA (I) =XW(I) 

89  (1)  =YW  { I ) 

12  CONTINUE  109 


*T  =  N 

CALL  PLOT SC(fl  At  ) 

FkiM  *  »**i  FLAG=  "t  TFL  AG 
50  U  11 

ccnt:  nlc 
'1 
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